## Abstract

We introduce a measurement (the *discrepancy*) of the minimum energy needed to transform from a standard parametrized planar curve *c*_{0} to an observed curve *c*_{1}. To this end, we say that a curve of transformations in the special Euclidean group SE(2) is *admissible* if it maps the source curve to the target curve under the point-wise action of SE(2) on the plane. After endowing the group SE(2) with a left-invariant metric, we define a relative geodesic in SE(2) to be a critical point of the energy functional associated to the metric, over all admissible curves. The discrepancy is then defined as the value of the energy of the minimizing relative geodesic. In the first part of the paper, we derive a scalar ODE which is a necessary condition for a curve in SE(2) to be a relative geodesic, and we discuss some of the properties of the discrepancy. In the second part of the paper, we consider discrete curves, and by means of a variational principle, we derive a system of discrete equations for the relative geodesics. We finish with several examples.

## 1. Introduction

Computational anatomy is the modern study of anatomical shape and its variability. This study originated in 1917, in the seminal book *On growth and form* by D'Arcy Thompson, who recognized that anatomical comparison is a mathematical problem, and that its solution would lie in what he called *the theory of transformations* [1, p. 1032].^{1}

Since then, D'Arcy Thompson has been proved correct, and many mathematical concepts, particularly concepts from Lie groups, Riemannian geometry and analysis have been applied to solve what is now called ‘the image registration problem’. An example is the comparison of medical images, which is a worldwide technology with an enormous number of uses every day. Perhaps not surprisingly, the comparison of medical images is still based on the theory of transformations, as enhanced by its modern developments.

Many mathematical frameworks have been developed to deal with the image registration problem, using the theory of transformations. The modern frameworks are well developed and fascinating for the geometry and analysis that underlies their solutions of the image registration problem. An outstanding example is the framework of *large deformation diffeomorphic metric mapping* (LDDMM).^{2} Many papers and books have been written about LDDMM and its variants. We refer the reader to the fundamental texts by Mumford & Desolneux [2] and Younes [3].

### (a) Contributions of this paper

This work arose in the context of LDDMM, but it reformulates the problem of registration in a simple and direct manner, for the comparison of two planar curves without using the diffeomorphism group. The reformulation is based on the idea that the *discrepancy* between two planar curves could be estimated quantitatively by defining the minimal amount of deformation needed to deform the source curve into the target curve using only the transformations of the special Euclidean group SE(2) of rotations and translations acting point-wise on the plane. The point of the paper is to define *discrepancy* precisely enough that the quantitative comparison will be meaningful.

Before explaining the main content of the paper and how it is organized, we present a few definitions and notation that set the context.

Let **c**_{0},**c**_{1}:[0,1]→*M* be curves, where *M* is a manifold. Let *G* be a Lie group with a Riemannian metric 〈 ,〉, and a transitive left action ⋅ on *M*. A function *g*: [0,1]→*G* is called *admissible* when *g*(*s*)⋅**c**_{0}(*s*)=**c**_{1}(*s*), for all *s*∈[0,1]. The *energy* of an admissible *g* is
where ∥ ∥_{g(s)} denotes the Riemannian norm. A critical point *g* of *E*_{0} is said to be a *geodesic relative to* **c**_{0},**c**_{1} or just a *relative geodesic* when **c**_{0},**c**_{1} are understood.

This paper investigates relative geodesics in the special case where *G* is the group SE(2) of Euclidean motions of , with a particular kind of left-invariant Riemannian metric, and the standard left action on . Then, given , their *discrepancy* *δ*(**c**_{0},**c**_{1}) is the minimum of *E*_{0}(*g*) as *g* varies over all admissible curves in SE(2). A minimizer *g* of *E*_{0} is necessarily a relative geodesic although, as seen in proposition 3.8, not all relative geodesics are minimizers.

When the discrepancy is 0, there is a constant admissible *g*, namely **c**_{0},**c**_{1} are congruent. More generally, the minimum energy *δ*(**c**_{0},**c**_{1}) is the minimum total variability of an admissible *g*. This is meant to capture, to some degree, the intuitive difficulty we find in transforming by eye from the parametric curve **c**_{0} to **c**_{1}. Note that the discrepancy depends on the parametrization of both curves, in contrast to LDDMM, so that we expect discrepancy to be the most useful in practice for curves which come with a natural parametrization, for instance samples of handwriting.

Because the Riemannian metric on SE(2) is chosen to be left-invariant, it follows easily that *δ*(**c**_{0},*h*.**c**_{1})=*δ*(**c**_{0},**c**_{1}) for any constant *h*∈*G*. Usually however *δ*(*h*.**c**_{0},**c**_{1})≠*δ*(**c**_{0},**c**_{1}), as in proposition 3.7. The discrepancy is therefore asymmetric in *c*_{0} and *c*_{1}, in the same way that, for most equal-length lists *x*_{0} and *x*_{1} of real numbers, the regressions of *x*_{1} on *x*_{0} and of *x*_{0} on *x*_{1} yield different linear relationships.

Sometimes there are continua of relative geodesics, as in lemmas 3.5 and 3.6. These facts are found by studying the Euler–Lagrange equations (called the continuous equations of motion) for relative geodesics, together with the natural boundary conditions at *s*=0,1. The Euler–Lagrange equations for relative geodesics are derived in §3 by Euler–Poincaré reduction, and also directly.

In §4, we introduce discrete analogues of relative geodesics in SE(2), as objects of separate interest, and with a view to developing numerical methods for continuous relative geodesics. Discrete curves are finite sequences, and the discrete energy is defined in terms of the Cayley map for SE(2). The Euler–Poincaré approach adapts nicely to the discrete case, leading to projected discrete equations of motion (4.22), and projected discrete boundary conditions (4.23), from which discrete relative geodesics can be found by Newton iteration. Alternatively, discrete relative geodesics can be calculated by direct minimization of the discrete energy. In §6, the numerical methods developed in §4 are used to illustrate properties of relative geodesics and discrepancies for some simple examples of plane curves. A morphing procedure, using a minimal relative geodesic, illustrates the geometrical difficulty of transforming **c**_{0} to **c**_{1}.

## 2. The Euclidean group SE(2)

In this section, we recall some basic results about the Lie group SE(2), which consists of rotations and translations in the Euclidean plane. For more information about SE(2) and its role in mechanics and control, see [4–6].

### (a) The Lie group SE(2)

The elements of SE(2) are pairs *g*:=(*R*_{θ},**x**), where *R*_{θ}∈*SO*(2) represents the counterclockwise rotation over an angle and represents the translation along **x**. There is a one-to-one correspondence between elements *g*=(*R*_{θ},**x**)∈SE(2) and 3×3 matrices of the form
2.1

In terms of matrices, the group multiplication in SE(2) is given by matrix multiplication. In components, we have that (*R*_{θ},**x**)⋅(*R*_{φ},**y**)=(*R*_{θ+φ},*R*_{θ}**x**+**y**) and (*R*_{θ},**x**)^{−1}=(*R*_{−θ},−*R*_{−θ}**x**).

### (b) The Lie algebra

The Lie algebra of SE(2) will be denoted by . The elements of the three-dimensional Lie algebra can be viewed as infinitesimal rotations and translations in the plane and are represented as column vectors , where and . There is a one-to-one correspondence between the elements *ξ* of and 3×3 matrices of the form
2.2In terms of matrices, the Lie bracket on is given by the matrix commutator: , for all . In components, we have that [(*ω*,**v**),(*η*,**w**)]=(0,−*ωJ***w**+*ηJ***v**).

The elements of the dual space are likewise column vectors in , denoted as *μ*=(*π*,**p**)^{T} with and . The duality pairing is given in terms of the Euclidean inner product on by
2.3

### (c) Choice of a norm on SE(2)

Let *m*>0 be a fixed parameter. We define a norm on the vector space by
2.4for .

We extend this norm by left translation to a norm on the tangent vectors to the group SE(2), given by
2.5for all *v*_{g}∈*T*_{g}SE(2). Here, we represent *g* in the matrix form (2.1) and the tangent vectors *v*_{g} as 3×3 matrices
where . The multiplication on the right-hand side of (2.5) is matrix multiplication, and . Upon expanding the matrix product, we find for the norm
2.6Note that the norm is by definition invariant with respect to the left action of SE(2) on itself,

From now on, we will drop the subscripts on the norms just defined, and we will denote both by ∥⋅∥.

### (d) The action of SE(2) on

The group SE(2) acts on the plane in the standard way: an element (*R*_{θ},**x**)∈SE(2) transforms a point into **c**′=*R*_{θ}**c**+**x**. This action translates into an infinitesimal action of on defined by (*ω*,**v**)⋅**c**=−*ωJ***c**+**v**. For fixed, we denote by the *isotropy subalgebra* of elements in that fix **c**, that is,

We let be the annihilator of in . In other words, consists of all covectors (*π*,**p**) which vanish when contracted with the elements of . A small calculation shows that
2.7Note that is isomorphic to , whereas is isomorphic to . Lastly, we define the projection by
2.8so that is precisely the kernel of . This map will be useful later on.

## 3. Continuous relative geodesics in SE(2)

Throughout this section, we let be two fixed parametrized curves. We say that a curve *g*:[0,1]→SE(2) is *admissible* with respect to **c**_{0} and **c**_{1} if
where the dot on the left-hand side represents the standard action of SE(2) on . In other words, a curve *g*(*s*)=(*R*_{θ(s)},**x**(*s*)) is admissible if
3.1

### (a) The deformation energy

We now wish to find the admissible curves in SE(2) that minimize the *deformation energy*,
3.2where the norm ∥⋅∥ was given in (2.5), and the prime ′ represents the derivative with respect to the curve parameter *s*. The deformation energy measures the change in *g*(*s*) as *s* varies.

### Definition 3.1

Let be two parametrized curves. An admissible curve *g*:[0,1]→ SE(2) with respect to **c**_{0},**c**_{1} is a *relative geodesic* if it is a minimum of the deformation energy *E*_{0} over all admissible curves.

As the norm (2.5) is invariant with respect to the multiplication from the left by elements of SE(2), we may write the deformation energy equivalently as
where in the second expression, the definition (2.4) has been used. Here, in the Lie algebra is related to the group element *g*(*s*)=(*R*_{θ(s)},**x**(*s*)) by means of the equations
3.3which follow from expanding the left-trivialized derivative *g*^{−1}(*s*)*g*′(*s*),
Using the identification (2.2), we then arrive at equation (3.3). These relations are referred to as the *reconstruction relations*: given *ω*(*s*) and **v**(*s*), (3.3) may be viewed as a set of first-order ODEs specifying the components of *g*(*s*).

For the purpose of deriving the equations that determine the extremals of *E*_{0}, it will be convenient to add the reconstruction relations as constraints to the deformation energy, so that we obtain
3.4Here, *E* depends now on the curve (*R*_{θ},**x**), the Lie algebra elements , and the Lagrange multipliers (*π*,**p**), which can be viewed as elements of the dual of the Lie algebra. It will be shown below that the critical points of this augmented functional coincide with the critical points of the original deformation energy, given in (3.2).

Note that *E* can be written in a more concise, Lie-algebraic way as
3.5where , , and the brackets 〈〈⋅,⋅〉〉 refer to the inner product associated to the norm (2.4). This energy function can be generalized in a straightforward way to the case of relative geodesics with values in an arbitrary Lie group *G*.

The variational principle for (3.5), in which the configuration variables, velocities and momenta are varied independently while the reconstruction equations are treated as constraints, is a particular example of the *Hamiltonian–Pontryagin principle*, as introduced by Yoshimura & Marsden [7]. A version of the Hamilton–Pontryagin principle specific to Lie groups can be found in Cendra *et al.* [8]; see also Bou-Rabee & Marsden [9]. Our variational principle is also related to the *Clebsch variational principle* of Cotter & Holm [10] and Gay-Balmaz & Ratiu [11], although it does not coincide with it. Other approaches to matching of shapes in Euclidean spaces may be found, for instance, in Srivastava *et al.* [12].

### (b) The continuous equations of motion

We now derive the differential equations that describe the critical points of the deformation energy. Because of the analogy with the Euler–Lagrange equations in mechanics, we will refer to these equations as equations of motion.

We take variations of the augmented deformation energy *E* in (3.4) with respect to the variables (*θ*,**x**,*ω*,**v**,*π*,**p**), where the velocities (*ω*,**v**) and the momenta (*π*,**p**) are varied freely, while the configuration variables (*θ*,**x**) are varied with respect to variations that preserve the admissibility constraint (3.1). On the level of the variations, the infinitesimal version of this constraint is given by
3.6where the matrix *J* was given in (2.2), and this relation allows us to eliminate the variation *δ***x** in terms of *δθ*. This infinitesimal constraint was obtained by taking a one-parameter family (*θ*_{ϵ}(*s*),**x**_{ϵ}(*s*)) of admissible curves in SE(2): taking the derivative of the admissibility constraint (3.1) with respect to *ϵ*, and putting
then gives (3.6).

Taking variations of *E* first with respect to *π* and **p** results in
so that if a curve *s*↦(*θ*(*s*),**x**(*s*),*ω*(*s*),**v**(*s*),*π*,**p**(*s*)) is a critical point of *E*, then the reconstruction equation (3.3) must hold. Taking variations with respect to the velocities *ω* and **v** similarly results in
so that for a critical point of *E*, the following Legendre transformations between the velocities and the momenta must hold:
3.7

Lastly, taking variations with respect to the configuration variables (*θ*,**x**), we obtain
where we have integrated by parts. We now use (3.6) to eliminate *δ***x** in the function of *δθ*, and we obtain
As *δθ* is arbitrary, we see that *δE* vanishes whenever the expressions preceding *δθ* on the right-hand side vanish, so that
3.8with boundary conditions *π*+**p**^{T}*J***c**_{0}=0 for *s*=0,1. Note that in these equations, *θ* and **x** are not arbitrary but are related by the admissibility constraint (3.1). By writing (3.8) in terms of the Lie algebra quantities (*ω*,**v**), and with some simple algebraic simplifications, we finally arrive at the following result.

### Theorem 3.2

*Let* *be two parametrized curves. An admissible curve (R*_{θ(s)}*,***x***(s)) in* SE(2) *is a critical point of the deformation energy E if and only if the following scalar equation holds:
*3.9*with natural boundary conditions mθ′+***v**^{T}*J***c**_{0}*=0 for s=0,1. Here,* **v** *is given in terms of θ and* **x** *by the reconstruction relation of (3.3), and the admissibility constraint (3.1) holds.*

As the relative geodesics are precisely the minima of the deformation energy *E*, restricted to the space of admissible curves, equation (3.9) is a necessary condition for a curve *g*:[0,1]→SE(2) to be a relative geodesic.

The equation of motion (3.9) can be written in a form that involves the angle *θ* only. To this end, we differentiate the admissibility constraint (3.1) with respect to *s* and multiply from the left by *R*_{−θ} to obtain **v**=*R*_{−θ}**c**_{1}′−**c**_{0}′+*J***c**_{0}*θ*′. One further differentiation leads to an expression for **v**′, and these two expressions can be substituted into (3.9) to give rise to the following non-autonomous second-order ODE for *θ*:
with boundary conditions for *s*=0,1.

Given the two curves **c**_{0} and **c**_{1}, the previous equations form a boundary-value problem for *θ*. Once *θ* is determined from these equations, the linear displacement **x** can be found from the admissibility constraint (3.1).

To conclude, we remark that the derivation of equation (3.9) is very similar to the derivation of the (constrained) Euler–Poincaré equations of classical mechanics (e.g. [4]), whereby a mechanical system on a Lie group *G* is described by a system of equations on the Lie algebra . Because of this analogy, we will refer to (3.9) as the Euler–Poincaré equations for relative geodesics.

#### (i) A direct derivation of the equations of motion

In §2b, we present an alternative derivation of equation (3.9), which does not use the Euler–Poincaré framework as in the previous section. This derivation is arguably somewhat more straightforward than the one presented earlier, and we will use the resulting Euler–Lagrange equations extensively in §3*c* below. The advantage of the Euler–Poincaré equations, however, is that they can easily be discretized, as we shall show in §4.

We begin by introducing the function . Note that is precisely the integrand of the deformation energy *E*_{0} in §3*a*. We now take the derivative of the admissibility constraint (3.1), and use the resulting equation to obtain an expression for **x**′. Upon substituting this expression into , we obtain a function ℓ that depends on *θ* and *θ*′, and is given by
3.10

The function ℓ can now be viewed as a Lagrangian function on the tangent bundle ; its Euler–Lagrange equations are 3.11

For further reference, we define the momentum conjugate to *θ* as
3.12and compute . By substituting these expressions into the Euler–Lagrange equation (3.11), we obtain yet another set of equations characterizing relative geodesics, which we summarize in the following theorem.

### Theorem 3.3

*Let* *be two parametrized curves. An admissible curve (R*_{θ(s)}*,***x***(s)) in* SE(2) *is a critical point of the deformation energy E*_{0} *if and only if the following Euler–Lagrange equation for θ holds:
*3.13*with natural boundary conditions given by p=0 for s=0,1. Here, p(s) is given by (3.12) and* **x***(s) is expressed as a function of θ(s) using the admissibility constraint (3.1).*

The procedure of substituting the constraints into the deformation energy to obtain a Lagrangian function that depends on fewer degrees of freedom is similar to the approach of Chaplygin for systems with non-holonomic kinematic constraints. In this approach, one eliminates the constrained degrees of freedom to obtain a system of reduced Euler–Lagrange equations with gyroscopic forces. The latter vanish if the constraints are integrable, as is the case for relative geodesics.

### (c) Discrepancy between planar curves

### Definition 3.4

Let be two parametrized curves in the plane. The *discrepancy* *δ*(**c**_{0},**c**_{1}) between **c**_{0} and **c**_{1} is the minimum of the deformation energy *E*_{0} over all (**c**_{0},**c**_{1})-admissible curves
where *Adm*(**c**_{0},**c**_{1}) is the set of all (**c**_{0},**c**_{1})-admissible curves.

Note that the admissible curve *g*:[0,1]→SE(2) which minimizes *E*_{0} can be found among the solutions of the equations of motion derived in §2b.

#### (i) Asymmetry of the discrepancy

The discrepancy *δ*(**c**_{0},**c**_{1}) provides a measure of the difference between the curves **c**_{0} and **c**_{1}. Here, we show that the discrepancy is in general not symmetric, that is, *δ*(**c**_{0},**c**_{1}) differs in general from *δ*(**c**_{1},**c**_{0}).

Throughout the remainder of §2c(i), we use the formulation of the Euler–Lagrange equations given in theorem 3.3.

### Lemma 3.5

*Let* **c**_{0}(*s*)=**0** *for all s. The geodesics relative to* (**c**_{0},**c**_{1}) *are parametrized by θ(0), and in each case*,

### Proof.

For a geodesic relative to (**c**_{0},**c**_{1}), we have that *p*=*mθ*′, and from the Euler–Lagrange equation (3.13), it follows that *mθ*′′=0, so that *θ* is an affine function of *s*. By using the natural boundary conditions, it follows that *θ* is constant. Conversely, any constant *θ* satisfies the Euler–Lagrange equation (3.13) and the natural boundary conditions, and therefore defines a relative geodesic. □

### Lemma 3.6

*Let* **c**_{1}(*s*)=**0** *for all s. The geodesics relative to* (**c**_{0},**c**_{1}) *are parametrized by θ*(0), *and in each case, the discrepancy is given by*
3.14

### Proof.

By (3.13), *p* is constant, and by the natural boundary conditions, *p* is identically 0, so that
By substituting this expression for *θ*′ into the deformation energy, we obtain (3.14) after some simplifications. □

### Proposition 3.7

*Let* **c**_{0}(*s*)=**0** *for all s. Then, δ*(**c**_{1},**c**_{0})≤*δ*(**c**_{0},**c**_{1}), *with equality only in the case where, for some* *and some* **c**_{1}(*s*)=*ϕ*(*s*)**x**_{0}.

### Proof.

Combining lemmas 3.5 and 3.6, we have that
so that *δ*(**c**_{1},**c**_{0})<*δ*(**c**_{0},**c**_{1}), except when (**c**_{1}′)^{T}*J***c**_{1} vanishes identically. In this case, **c**_{1}′(*s*)=*μ*(*s*)**c**_{1}(*s*) for some function *μ*, and therefore **c**_{1}(*s*)=*ϕ*(*s*)**x**_{0}, with *ϕ*(*s*)=*e*^{μ(s)}. □

As an illustration, we take **c**_{0}(*s*)=**0** and for *s*∈[0,1]. By lemmas 3.5 and 3.6, we have that *δ*(**c**_{0},**c**_{1})=*π*^{2}/2≈4.93, and , so that indeed *δ*(**c**_{0},**c**_{1})>*δ*(**c**_{1},**c**_{0}).

#### (ii) Non-minimizing relative geodesics

There are always at least two geodesics relative to (**c**_{0},**c**_{1}), corresponding to critical points of the deformation energy *E*_{0} regarded as a function of *θ*_{0} where the relative geodesic *s*↦(*R*_{θ(s)},**x**(*s*)) satisfies (3.13) for all *s* and the natural boundary conditions at *s*=0.

### Proposition 3.8

*Let* **c**_{1} *be a non-constant affine line segment, and suppose* **c**_{0}(0)=**0** *with* **c**_{0}(1)≠**0**. *Then there are exactly two geodesics relative to* (**c**_{0},**c**_{1}), *and these are determined by*
*with χ*_{1} *the angle between* **c**_{0}(1) *and* **c**_{1}′(1). *Only one of these relative geodesics is a global minimizer of E*_{0}.

### Proof.

As **c**_{1}′(*s*)=:**C**, where **C** is a constant vector, the right-hand side of the Euler–Lagrange equation (3.13) can be written as a total *s*-derivative,
Consequently, the Euler–Lagrange equations imply that the quantity
3.15is conserved. As **c**_{0}(0)=**0** and *p*=0 at *s*=0, is identically 0. So, for some *θ*_{0},
At the terminal end of the curve, i.e. for *s*=1, we have that . From (3.12), it then follows that , namely *R*_{θ(1)}**c**_{0}(1) is a multiple of **C**, so that *θ*(1)=±*χ*_{1} for some angle *χ*_{1}. As a consequence, the initial angle *θ*_{0} satisfies
Considering *E*_{0} as a function of , one of these values of *θ*_{0} is a point of global maximum, the other is a point of global minimum, and *E*_{0} has no other critical points. □

As an illustration, we consider the discrepancy between two line segments. We take **c**_{0}(*s*)=*s***e**_{x} and **c**_{1}(*s*)=*s***e**_{y}, with **e**_{x},**e**_{y} the standard unit vectors along the positive *x*- and *y*-axes, respectively. From the previous proposition, we deduce that *θ*(*s*)=±*π*/2 for all *s*. For the solution with *θ*(*s*)=*π*/2, the admissibility condition results in **x**(*s*)=0 for all *s*. In this case, the effect of applying the relative geodesic is to rotate all of the points of **c**_{0} over *π*/2 and not effect any translation. As the relative geodesic is constant, *g*(*s*)=(*R*_{π/2},**0**), the deformation energy vanishes identically, so that *g*(*s*) is a minimizing geodesic. In the case where *θ*(*s*)=−*π*/2, the admissibility constraint yields **x**(*s*)=2*s***e**_{y} so that the effect of the relative geodesic is to rotate each point **c**_{0}(*s*) over −*π*/2, followed by a translation over 2*s***e**_{y}. The deformation energy in this case is *E*=2.

### Remark 3.9

In the proof of proposition 3.8, we have seen that the quantity is conserved when **c**_{1}′(*s*)=**C** with **C** constant. A natural question to ask is the following: is there a continuous symmetry whose associated conserved quantity (through Noether's theorem) is precisely ? To see that this is indeed the case, we return to the Lagrangian ℓ in (3.10), which we rewrite as
where is defined as
3.16

As ℓ and differ by a total *s*-derivative, they give rise to the same Euler–Lagrange equations [13]. For , the momentum conjugate to *θ* is precisely the quantity defined in (3.15),
In the case that **c**_{1}′′=0, we see from (3.16) that does not depend on *θ*, so that is a conserved quantity,

## 4. Discrete relative geodesics in SE(2)

We now assume that we have two discrete curves (**c**_{0})_{k},(**c**_{1})_{k}, *k*=0,…,*N* of *N* points each. We wish to find a discrete curve *g*_{k}=(*R*_{θk},**x**_{k}), *k*=0,…,*N* in SE(2) which is *admissible* in the sense that
4.1for all *k*=0,…,*N*. To derive a discrete version of the deformation energy *E*, we need to discretize the spatial derivatives that appear in (3.5). We do this by means of the Cayley map from to SE(2).

Our way of discretizing the variational principle, as well as the discrete equations obtained from it, is inspired by the *discrete Hamilton–Pontryagin principle* of Bou-Rabee & Marsden [9]; see also Kobilarov [14] and Stern [15]. As in the continuous case, the main difficulty here is the incorporation of the admissibility constraint (4.1).

### (a) The Cayley map

### Definition 4.1

The Cayley map is given by
4.2where, if ,
4.3with ** v**=(

*v*

_{1},

*v*

_{2}). Note that depends only on

*ω*and is in fact the Cayley transform in

*SO*(2).

The Cayley map is in fact a (1,1)-Padé approximation to the exponential map from . In contrast to the exponential, the Cayley map has the advantage that it is an algebraic map, so that it is easily computable. The Cayley map shares with the exponential map a number of useful properties, which will be used in some of the derivations below 4.4for all , where id. is the identity matrix.

#### (i) The right-trivialized derivative

For our purposes, we will need the *right-trivialized derivative* of the Cayley map, defined by Iserles *et al.* [16],
4.5Note that *D*Cay(*ξ*)⋅*η* is an element of *T*_{ξ}SE(2), which is translated back to by right-multiplication by Cay(*ξ*)^{−1}. In this way, dCay is a map from to , which is linear in the second argument.

For fixed , we denote the inverse of dCay_{ξ} by . From (4.5), we have that
4.6

For the group SE(2), the Cayley map and its derivatives were computed explicitly in Kobilarov [14]. Keeping in mind that the elements of are represented as column vectors, for each , is a linear transformation from to itself, given by , where, for *ξ*=(*ω*,*v*,*w*), the matrix *M*(*ξ*) is given by
4.7

For future reference, we record the following property of the right-trivialized derivative [9], for all :
4.8where Ad is the adjoint action of SE(2) on , defined by Ad_{g}(*ξ*)=*gξg*^{−1}, where the elements on the right-hand side are interpreted as matrices, as in (2.1) and (2.2).

Lastly, for each , its adjoint, , is a linear map from to itself, defined by 4.9relative to the duality pairing (2.3). Explicitly, 4.10

We will use the Cayley map to provide a parametrization of a neighbourhood of the identity in SE(2) by means of the Lie algebra , but it is possible to replace the Cayley map by any other local diffeomorphism satisfying (4.4) from to SE(2), such as the exponential map. The Cayley map, however, has the advantage that it is efficiently computable, and its derivative is particularly easy to characterize.

### (b) The deformation energy

#### (i) The discrete reconstruction relations

Using the Cayley map, we discretize the reconstruction relations (3.3) as follows. Given two successive elements *g*_{k},*g*_{k+1} in SE(2), we define the *update element* by
4.11This is the discrete counterpart of the relation *ξ*=*g*^{−1}*g*′. Explicitly, if *ξ*_{k}=(*ω*_{k},**v**_{k}) and *g*_{i}=(*R*_{θi},**x**_{i}), *i*=*k*,*k*+1, we have for the components
4.12The first relation is equivalent to the following trigonometric relation:
4.13

#### (ii) The deformation energy

To discretize the deformation energy, we now proceed as in the continuous case. We define *E* as
4.14which can be viewed as the discrete counterpart of (3.5). Here, *g*_{k}∈SE(2), and are independent variables. As mentioned at the beginning of §4, this energy function was originally introduced in Bou-Rabee & Marsden [9], and the derivations up to (4.17), when we have to enforce the discrete admissibility constraint, will follow that paper.

By taking variations with respect to *μ*_{k}, we recover the definition (4.11) of the update element *ξ*_{k}. By taking variations with respect to *ξ*_{k}, we arrive at the equation , or explicitly
4.15

Lastly, by taking variations with respect to the group element *g*_{k}, we obtain
where we have used the fact that . We now introduce the quantity and focus first on the first derivative term, which we write as
where we have used the definition (4.6) of dCay^{−1}, together with the expression (4.11) for the update element. For the second term in *δE*, we proceed along similar lines,
where we have used the property (4.8) in the last step.

Substituting both of these expressions for the derivatives back into the expression for *δE*, we arrive at
where we have introduced the adjoint of the linear map ( using the definition (4.9). We now rearrange the terms in the sum to get
4.16

It remains for us to obtain an expression for the variations . As *E* is varied over all discrete admissible curves, (4.1) must hold, and by differentiating and multiplying by *R*_{−θk}, we find
4.17where *δθ*_{k} and **w**_{k} are the components of *σ*_{k}. Note that if *δg*_{k}=(*δθ*_{k},*δ***x**_{k}), then **w**_{k}=*R*_{−θk}*δ***x**_{k}. In other words, *σ*_{k} is an element of . As *σ*_{k} is otherwise arbitrary, we arrive at the following weak form of the discrete equations of motion:
4.18for all , together with the weak boundary conditions
4.19for all and . Another way to formulate the equations of motion is to observe that the left-most factor in each of these contractions must take values in , defined in (2.7). In this way, we arrive at the following theorem.

### Theorem 4.2

*A discrete admissible curve g*_{k}*∈SE(2), k=0,…,N, is a critical point of the deformation energy E if and only if (4.18) holds, with boundary conditions (4.19). This is equivalent to
*4.20*together with the boundary conditions
**Here, the discrete reconstruction equation (4.12) holds, and the discrete momenta are given by* *as in (4.15).*

#### (iii) The discrete constraint

Equation (4.20) is a single scalar equation, which in itself is insufficient to determine all three components of *ξ*_{k}=(*ω*_{k},**v**_{k}). We now show that the admissibility condition (4.1) gives rise to two further equations, allowing all three components of *ξ*_{k} to be determined.

By taking the admissibility constraint for *k* and subtracting it from the constraint for *k*+1, we obtain (after multiplying from the left by *R*_{−θk}) that
Using the relation (4.12) for the components of the Cayley map, this becomes
4.21Given *ω*_{k}, the first component of *ξ*_{k}, this relation can be solved to find the corresponding linear velocity **v**_{k}. In fact, as depends linearly on **v**_{k}, (4.21) is just a linear equation for **v**_{k}.

*Summary.* To convince ourselves that the equations derived so far are sufficient to determine the discrete curve *g*_{k}, *k*=0,…,*N* completely, we summarize the equations of motion. A practical way to implement these equations will be given in §4*c*.

Assume that two successive elements *g*_{k−1},*g*_{k}∈SE(2) are given, which satisfy the admissibility condition (4.1). The equations allow for *g*_{k+1} to be found as follows:

— using the Cayley transform (4.11), find ;

— solve the following two equations simultaneously for

*ξ*_{k}=(*ω*_{k},**v**_{k}): (4.20) and (4.21); and— given

*ξ*_{k}and*g*_{k}, determine*g*_{k+1}from the update relation (4.11).

### (c) Practical implementation of the discrete equations of motion

Using the projector defined in (2.8), we first write the equations of motion as
4.22and the boundary conditions as
4.23Here, we use the matrix expression *M*(*ξ*) for dCay^{−1}, given by (4.7) and , . From now on, we will drop the hat over the linear quantities *ξ* and *μ*, as the factors of *h* in front of *ξ* and *μ* can be restored at a later stage without any ambiguity.

#### (i) The discrete equations of motion for SE(2)

We introduce the matrices
4.24and observe that
Using the fact that *B*(*ω*)*J*=*J*−(*ω*/2)*I*, we obtain for the projection (2.8) that for any point **c**,
4.25

Using (4.15) to write the momenta in terms of *ω* and **v**, we obtain, for the projected equation of motion (4.22),
4.26

#### (ii) The linear velocities

To solve the discrete constraint (4.21) for **v**, we let
4.27so that (4.21) is equivalent to *A*(*ω*_{k})**v**_{k}=**b**_{k}, which can be solved as
4.28Notice that the right-hand side depends only on *ω*_{k} and *θ*_{k}.

#### (iii) The boundary conditions

Lastly, we show how the boundary conditions (4.23) can be made more explicit. We assume that the two discrete curves have been translated to the origin, so that (**c**_{0})_{0}=(**c**_{1})_{0}=0. Using (4.25), the boundary condition for *k*=0 becomes , so that *ω*_{0}=0.

To find **v**_{0}, we focus on the discrete constraint (4.21) for *k*=0, which becomes , so that at *k*=0, the following conditions hold:
4.29where *θ*_{0} is arbitrary. Once *θ*_{0} is chosen, these relations suffice to find the first two group elements *g*_{0} and *g*_{1}. At the other end of the curve, the boundary condition (4.23) for *k*=*N* reads
4.30

### (d) Solving the boundary-value problem

We now outline a procedure for solving the boundary-value problem (4.22), (4.23) by using a simple shooting method based on Newton iteration.

### (e) First-variation equations

We begin by linearizing equation (4.26) around a given solution. We denote
4.31and **d**_{k}(*ω*):=(*ω*/2)**v**+*J*(**c**_{0})_{k}−(*ω*/2)(**c**_{0})_{k}. The first-variation equation may then be expressed as
4.32This is a single linear equation for *δω*_{k}; *δ***v**_{k} can be obtained from the linearization of (4.28). After some algebra, we obtain
where *A*(*ω*) was defined in (4.24). Lastly, the variation *δθ*_{k} may be obtained from the linearization of the Cayley equation (4.13) and is given by
The initial conditions for the variational equations can be obtained by linearizing (4.29),
4.33while *δθ*_{0} is arbitrary.

### (f) Newton iteration

The shooting method now proceeds as follows. Starting with a value for the angle *θ*_{0}, we solve the initial value constraints (4.29) for (*ω*_{0},**v**_{0}) and the equations of motion (4.26) to obtain a discrete solution trajectory (*ω*_{k},**v**_{k}) for *k*=1,…,*N*−1. At the terminal end, the boundary conditions (4.30) will in general not be satisfied, so that we need to adjust the initial angle *θ*_{0}. We do this by means of a Newton iteration.

We denote by *Φ*(*θ*_{0}) the value of the terminal boundary condition (4.30), evaluated on the solution trajectory corresponding to *θ*_{0}. If *Φ*(*θ*_{0})≠0, that is, *θ*_{0} leads to a solution that does not satisfy the terminal boundary condition, an improved value for *θ*_{0} is given by *θ*_{0}↦*θ*_{0}−*Δθ*_{0}, where
The derivative ∂*Φ*/∂*θ* of the terminal boundary condition with respect to the initial condition can be computed as follows:
where is the right-hand side of (4.30), viewed as a function of *ω*_{N−1} and **v**_{N−1}, and where the derivatives of have been expressed in terms of *f*_{N} and **d**_{N} defined in (4.31). The derivatives ∂*ω*_{N−1}/∂*θ* and ∂**v***ω*_{N−1}/∂*θ* of the terminal quantities with respect to the initial condition *θ* can be expressed in terms of the variational equations: if we solve the variational equations with initial conditions *δθ*_{0}=1, *δω*_{0}=0 and *δ***v**_{0} given by (4.33), then
Putting all of this together, we then arrive at the following expressing for *Δθ*_{0}:
4.34

If we view the energy *E* of a solution as a function of the initial angle *θ*_{0}, then we are effectively searching for the minimum of a scalar function defined on a compact set. The minimum is therefore guaranteed to exist and is typically found in practice after only a few iterations of the algorithm.

## 5. Direct minimization of the energy functional

Instead of explicitly solving the boundary-value problem (4.22), one can also minimize the deformation energy (4.14) directly. To bring *E* into a form that can be handled conveniently by standard optimization software, we write it as
5.1Here, we recall that the linear quantities *ω*_{k} and **v**_{k} have been scaled by *h*; this explains the factor *h* in front of *E* on the left-hand side. Recall that *ω*_{k} is given in terms of the angles *θ*_{l} by (4.13), whereas **v**_{k} is given by (4.28). Equivalently, *E*_{0} can be written as
where **b**_{k} is given by (4.27).

The gradient of *E*_{0} with respect to *θ*_{i} can easily be computed from this expression. A standard calculation yields
for *k*=1,…,*N*−1, where
and , while **g**_{i}=*JR*_{−θi}((**c**_{1})_{i+1}−(**c**_{1})_{i}). At the terminal points, we have and ∂*hE*_{0}/∂*θ*_{N}=*Q*_{N−1}.

## 6. Numerical results

We have implemented the algorithm of §5 in Python programming language. Our implementation, together with scripts for the simulations in this section, can be found at https://github.com/jvkersch/se2-curve-matching and in the electronic supplementary material.

Unless otherwise noted, the parameter *m* in the deformation energy (4.14) will be taken to be equal to 2 in the numerical experiments below.

### (a) Discrepancy of simple planar curves

In this example, we compute the discrepancy between the unit circle with parametric representation and the figure eight given by , for *s*∈[0,1]. We divide the parameter interval into *N*=100 equal subintervals, and we let (**c**_{i})_{k}=**c**_{i}(*hk*) for *i*=0,1 and *k*=0,…,*N*, where *h*=1/*N*. We finish by rotating and translating both curves so that they start at the origin and are tangential to the *x*-axis. The final layout of the curves is shown in figure 1. Upon solving the boundary-value problem, we find that the global minimum of the deformation energy is located at , where .

Last, we use the problem of finding the discrepancy between the circle and the figure eight to get a rough estimate of the order of accuracy of our numerical method in terms of the number of sample points on the discrete curves. To do this, we run a number of simulations: at each stage, we choose *N*+1 sample points on each curve, and we compute the discrepancy between **c**_{0} and **c**_{1}. As *N* increases, the discrepancy is expected to approach a limit value, and the rate at which convergence takes place will give us an estimate of the order of accuracy of our numerical method.

In table 1, we have listed the discrepancy for a few choices of *N*. Roughly speaking, as *N* goes up by a factor of 10, the discrepancy gains two digits of accuracy, so the order of accuracy of the method is approximately 1.

### (b) Interpolation for curves on SE(2)

To provide some further visual insight into the nature of relative geodesics, we use a simple form of interpolation on SE(2). Given an element *g*∈SE(2), we define, for *ϵ*∈[0,1],
where is the exponential map and its inverse. The element *g*_{ϵ} will be well defined provided that *g* is in the range of the exponential map. In this way, we obtain a smooth curve in SE(2) that connects the identity element to *g* as we let *ϵ* range from 0 to 1.

Given a curve *s*↦*g*(*s*), we may apply this transformation to each point of the curve in order to obtain a family of curves *g*_{ϵ}(*s*). Now, assume that the original curve *g*(*s*) matches the planar curves **c**_{0} and **c**_{1} and set **c**_{ϵ}(*s*)=*g*_{ϵ}(*s*)⋅**c**_{0}. As *ϵ* ranges from 0 to 1, **c**_{ϵ} will be a curve that smoothly interpolates between **c**_{0} and **c**_{1}. A similar procedure may be done for the case of discrete curves.

Most of the computations in the remainder of §6 were visualized using the following form of linear interpolation: whenever we compute an admissible curve *g*(*s*) with respect to two planar curves **c**_{0} and **c**_{1}, we construct the intermediate curves **c**_{ϵ} to give an idea of the deformations carried out by *g*(*s*). While the sequence of curves thus obtained has no immediate physical meaning, it nevertheless gives a good intuitive idea of the amount of deformation needed to match one curve with another. To illustrate this, we show in figure 2 the matching between a circle and a figure-eight shape, where the matching is first done using the global minimizer of the energy (3.2), and secondly with two different local minimizers.

### (c) Asymmetry of the discrepancy

We have mentioned before that the discrepancy is not symmetric. We illustrate this by matching a circle **c**_{0} of radius *r*=0.1 with a figure-eight shape, given by the parametric representation . On both curves, 100 points were sampled equidistantly in *s*, and the parameter *m* in the norm (3.2) was set to 2. We find the optimal admissible curve *g* by means of the algorithm in §4*d*. For the discrepancy, we obtain *δ*(**c**_{0},**c**_{1})=47.6261 and *δ*(**c**_{1},**c**_{0})=39.8011.

The deformations leading to these respective discrepancies are visualized in figure 3, using the interpolation procedure described in §6*b*.

As these deformations appear quite similar, despite the large difference in the discrepancies, we investigate some further characteristics of the minimizer. In figure 4, we plot the *θ*-coordinate of the minimizer for both matching problems. These curves appear quite different.

### (d) Discrepancy of polynomial curves

In this last example, we take a closer look at the relation between discrepancy and other geometrical invariants, in particular the total absolute curvature *κ*, defined as the integral of the absolute value of the curvature: .

We focus on unit-length segments of the polynomial curves *y*=*x*^{p}, for *p*=1,2,…, which we denote by **c**_{p}(*s*)=(*x*_{p}(*s*),*x*_{p}(*s*)^{p}), where the parameter *s*∈[0,1] denotes arclength with *s*=0 corresponding to the origin, and so the curve is traversed in the direction of the positive *x*-axis. A few of these curves are plotted in figure 5. For each of these curves **c**_{p}(*s*), we choose *N* equidistant points *s*_{k} in the parameter interval [0,1], so that *s*_{k}=(*k*−1)/(*N*−1) for *k*=1,…,*N*, and we set (**c**_{p})_{k}=**c**_{p}(*s*_{k}). In this way, we obtain for each exponent *p* a discrete curve consisting of *N* points (**c**_{p})_{k}, for *k*=1,…,*N*, at equal distance (along the curve) from each other.

For each exponent *p*>0, we compute the discrepancy between the discrete curve (**c**_{p})_{k} and the fixed curve (**c**_{0})_{k}, which is parallel to the *x*-axis. We have tabulated the results for a selection of exponents *p* in table 2, along with the total geodesic curvatures for each of the curves **c**_{p}. In figure 5, both invariants have been plotted for *p*=1,2,…,25. We see that the geodesic curvature increases as a function of *p*, corresponding to the more pronounced bend in the curve for higher *p*. On the other hand, the discrepancy first increases until *p*=6 and then decreases: while the curves for high *p* are more curved, the curvature is more localized and the curves as a whole are close to the *x*-axis, resulting in lower discrepancy.

## 7. Conclusions and outlook

In this paper, we have outlined a new measure for the discrepancy between planar curves, based on deforming one curve into the other by means of parameter-dependent transformations with values in the Lie group SE(2). We defined a relative geodesic in SE(2) to be a curve of transformations that extremizes a certain energy functional while mapping the first curve into the second, and we defined the discrepancy to be the value of the energy associated to the minimizing relative geodesic.

One of the advantages of our approach is that it can be generalized in a straightforward manner to deal with, for instance, discrepancies and relative geodesics between other types of geometric objects, such as curves in three-dimensional or two-dimensional images. Another direction for future research addresses the choice of Lie group of transformations. In this paper, we considered the group SE(2) of rotations and translations, but other groups acting on the plane can be treated similarly. For instance, one could imagine acting on the curves by means of shearing transformations and translations, and in this case, the relevant group is the semi-direct product , where SL(2) is the group of 2×2 matrices with unit determinant.

## Funding statement

D.D.H. and J.V. gratefully acknowledge partial support by the European Research Council Advanced grant no. 267382 FCCA. L.N. is grateful for support from this grant during a visit to Imperial College. J.V. is grateful for partial support by the IRSES project GEOMECH (no. 246981) within the 7th European Community Framework Programme and is on leave from a Postdoctoral Fellowship of the Research Foundation Flanders (FWO-Vlaanderen).

## Acknowledgements

We would like to thank Jaap Eldering, Henry Jacobs and David Meier for stimulating discussions and helpful remarks.

## Footnotes

- Received May 8, 2013.
- Accepted June 28, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.