M\"obius invariants of shapes and images

Identifying when different images are of the same object despite changes caused by imaging technologies, or processes such as growth, has many applications in fields such as computer vision and biological image analysis. One approach to this problem is to identify the group of possible transformations of the object and to find invariants to the action of that group, meaning that the object has the same values of the invariants despite the action of the group. In this paper we study the invariants of planar shapes and images under the M\"obius group $\mathrm{PSL}(2,\mathbb{C})$, which arises in the conformal camera model of vision and may also correspond to neurological aspects of vision, such as grouping of lines and circles. We survey the computational requirements of an invariant, and the known M\"obius invariants, and then develop an algorithm by which shapes can be recognised that is M\"obius- and parametrization-invariant, numerically stable, and robust to noise. We demonstrate the efficacy of this new invariant approach on sets of curves, and then develop a M\"obius-invariant signature of grey-scale images.

to its importance on fundamental grounds-it is one of the very few classes of Lie groups that act on the plane, it crops up in numerous branches of geometry and analysis, it is the smallest nonlinear planar group that contains the direct similarities, and it is the set of biholomorphic maps of the Riemann sphere-it also has direct applications in image processing since it arises in the conformal camera model of vision, in which scenes are projected radially onto a sphere [21,38,39]. It may also correspond to neurological aspects of vision, such as grouping of lines and circles (which are equivalent under Möbius transformations) [40].
From an applications point of view, different objects may be related by Möbius transformations, as is explored by Petukhov [33] for biological objects in fascinating detail in the context of Klein's Erlangen program 1 , giving examples of many body parts that are loxodromic to high accuracy (loxodromes have constant Möbius curvature and play the role in Möbius geometry that circles do in Euclidean geometry), that change shape by Möbius and other Lie group actions, and that grow via Möbius transformations (including 2D representations of the human skull). Other examples of 1-dimensional growth patterns, such as antenatal and postnatal human growth, appear to be well modelled by 1-dimensional linear fractional transformations x → (ax+b)/(cx+d). Discussing this work, Milnor [26] argued that "The geometrically simplest way to change the relative size of different body parts would be by a conformal transformation. It seems plausible that this simplest solution will often be the most efficient, so that natural selection tend to choose it." (Milnor was thinking of 3D conformal transformations, whose restriction to 2D is the Möbius transformations.) In this paper we present an integral invariant for the 2D Möbius group that is reparameterization independent, and demonstrate its use to identify curves that are related by Möbius transformations. We then consider the case of images, and describe an invariant signature by which Möbius transformed images can be recognised.
We begin in Section 2 by discussing the desirable properties of invariants for object classification, and introduce the property of bounded distortion, which subsumes most of the requirements identified in the literature. This is followed by an overview of the methods that have been used to give object classification invariants for curves and images (most often with respect to the Euclidean, similarity, and affine groups).
In Section 3 we present the classical Möbius invariants and discuss their utility with respect to the properties identified in Section 2, before using a numerical example based on an ellipse to demonstrate the difference between them. This is followed by the introduction of the invariant that we have identified as the best behaved for curves with respect to the requirements previously discussed. In order to demonstrate its utility, we present an experiment where a set of smooth Jordan curves are created, and then the invariant distance is computed, and compared with direct registration of each pair of curves in the Möbius group using the H 1 similarity metric. The results show that the invariant is well-behaved with respect to noise, and can be used to separate all but the most similar shapes.
In Section 4 we move on to images and demonstrate the use of a 3D Möbius signature that is very sensitive and relatively cheap to compute, while still being more robust than the analogous 1 D'Arcy Wentworth Thompson [37] famously deformed images of one species to match those of another; his theory of transformations is reviewed and interpreted in light of modern biology in [5]. In particular, it has given rise to image processing techniques such as Large Deformation Diffeomorphic Metric Mapping (LDDMM) [16], in which images are compared modulo infinite-dimensional groups such as the diffeomorphism group. Petukhov writes that Thompson "did not use the Erlanger program as the basis in this comparative analysis." However, the totality of Thompson's examples and the explanations in his text do indicate that in all cases he selected his transformations from the simplest group that would do an (in his view) acceptable job. Eight different groups are identified in [23, Table 1]; four are finite dimensional. Thus, whether consciously or not, Thompson's work was fully consistent with the Erlangen program. Of relevance to the present paper is that many of his examples use conformal mappings, and thus may be approximated (or even determined by) Möbius transformations. signature for curves. As far as we know, such differential invariant signatures (in the sense of Olver [30]) for images have not been considered previously.
2 Invariants of shapes and images 2

.1 Computational requirements of invariants for object classification
The mathematical definition of an invariant, namely, a G-invariant function on M , is not sufficiently strong for many computational applications. For example, object classification via invariants involves comparing the values of the invariant on different G-orbits, and the definition says nothing about this. Recognising this, Ghorbel [15] has given the following partial list of qualities needed for object recognition: (i) fast computation; (ii) good numerical approximation; (iii) powerful discrimination (if two objects are far apart modulo G, their invariants should be far apart); (iv) completeness (two objects should have the same invariants iff they are the same modulo G); (v) provide a G-invariant distance on M ; and (vi) stability (if two objects have nearby invariants, they should be nearby modulo G).
Calabi et al. [10] include in (ii) the requirement that the numerical approximation itself should be G-invariant, while Abu-Mostafa et al. [1] add the further desirable qualities of: (vii) robustness (if two objects are nearby modulo G, especially when one is a noisy version of the other, their invariants should be nearby); (viii) lack of redundancy (i.e. all invariants are independent); and (ix) lack of suppression (in which the invariants are insensitive to some features of the objects).
To this already rather demanding list we add one more, that the set of invariants should be: (xi) small.
The motivation for this last is purely a parsimony argument; it is easy to produce large sets of invariants without adding much utility, and smaller sets of easily computable invariants are to be preferred. This is particularly the case since some of these criteria are in conflict, so there will usually not be one invariant that is preferred for all applications; instead, the best choice will depend on the dataset and the particular application. Others are closely related, especially discrimination, distance, stability, robustness, and suppression; these all depend on which features of the objects are deemed to be signal and which are deemed to be noise.
These criteria can be unified and quantified in the following way: suppose that a distance on objects has been chosen that measures the features that we are interested in (the signal), and  Figure 1: Sketch of the approximate relationships between d G (x, y), which measures the distance between two objects x and y modulo a symmetry group G, and I(x) − I(y) , which measures the distance between their invariants, illustrating the role of the desired properties of completeness, robustness, stability, and discrimination. In object recognition, 'good discrimination' is sometimes taken to mean the avoidance of type I errors (relative to the null hypothesis that two objects are in the same group orbit, so that in a type I error dissimilar objects are classified as similar), the avoidance of type II errors (similar objects being classified as dissimilar), or is not specified; we take it to mean that both type I and type II errors are avoided. that is small for differences we are not interested in (the noise). This distance induces a distance on objects modulo G by (where · d is some appropriately chosen norm): (In principle one can compute d G (x, y) by optimizing over G, and indeed this is done in many applications. However, optimization is relatively difficult and unreliable, and becomes impractical on large sets of objects; this is the step that invariants are intended to avoid.) Suppose that a norm on invariants I : M → R k has been chosen. Then we can express the criteria as follows: These may be subsumed and strengthened by the property of bounded distortion: an invariant has bounded distortion with respect to d G if there exist positive constants c 1 , c 2 such that for all x, y ∈ M we have: Suppose the invariant is to be used to test the hypothesis that two objects are the same modulo symmetry group G. The situation is illustrated in Figure 1, which plots the distance between the invariants of two objects against the distance between the objects modulo G. The smaller the value of c 2 , the more robust the invariant is, and the fewer false positives will be reported. The larger the value of c 1 , the more stable the invariant is, the better its discrimination, and the fewer false negatives will be reported. If (1) holds only with c 1 = 0, the invariant is not complete; if there is no c 2 such that (1) holds, the invariant is not stable. In addition, in some applications one may be more interested in some parts of the space shown in Figure 1 than others: for example, in discriminating very similar or very dissimilar objects.
In practice, it may not be practical to find invariants with bounded distortion. Even completeness, the centrepiece of the mathematical theory of invariants, can be very demanding. It turns out that many invariants that are roughly described as 'complete' are not fully complete, that is, they do not distinguish all group orbits, but only distinguish almost all group orbits. For example, Ghorbel's Euclidean invariants of grey-level images [15] only distinguish images in which two particular Fourier coefficients are non-zero. Thus it will be stable only on images in which these two Fourier coefficients are bounded away from zero. The most commonly used Euclidean signature of curves, (κ, κ s ), is complete only on nondegenerate curves [18]. The fundamental issue is that even for very simple group actions, the space of orbits can be very complicated topologically, and thus very difficult to coordinatize via invariants. We illustrate these ideas via an example: Example 1 Consider n points z 1 , . . . , z n in the complex plane. Let G = S 1 act by rotating the points, ie. e iθ · (z 1 , . . . , z n ) = (e iθ z 1 , . . . , e iθ z n ). Then {z i z j : 1 ≤ i ≤ j ≤ n} forms a complete set of invariants. It is very large (there are n 2 real components) compared to the dimension 2n − 1 of the space, C n /S 1 , in which we are trying to work. However, if any one real component of the set is omitted, the resulting set is not complete. (For example, if |z 1 | 2 is omitted, then the points (1, 0, . . . , 0) and (2, 0, . . . , 0) have the same invariants, but do not lie on the same orbit.) This is the situation considered in the problem of phase retrieval [7]. By choosing combinations of {z i z j } it is possible to create smaller sets of complete invariants, but the computational complexity of calculating the description of the n points is still O(n 2 ) [7].
In cases where one settles for an incomplete set of invariants, or a set that is not of bounded distortion-so that the worst-case behaviour of the invariants is arbitrarily poor-it can make sense instead to study the average-case behaviour of the invariants over some distribution of the objects. This is analogous to the role that ill-posed and ill-conditioned problems play in numerical linear algebra, although in that case the ill-conditioning is intrinsic rather than imposed by considerations of computational complexity. Indeed, one can sometimes very usefully describe objects based on an extremely small number of invariants, for example, describing planar curves by their length and enclosed area, or by the first few Fourier moments of their curvature. What is wanted is to optimize the behaviour of the invariants, over some distribution of the objects, with respect to their time and/or space complexity. However, we know of no genuine cases in which such a program has been carried out.

Invariants of curves
There is a large literature concerning invariants to the Euclidean and similarity groups for both curves and images. The purpose of this section is to provide an overview of the dominant themes within that work that are relevant to our goal of identifying invariants for the Möbius group that satisfy at least some of the criteria listed in the previous section.
We define a (closed) shape as the image of a function φ : S 1 → R 2 . For different restrictions on φ this defines different spaces: if φ is a continuous injective mapping then this defines simple closed curves, while if φ is differentiable and φ (t) = 0 for all t this defines the 'shape space' Imm(S 1 , R 2 )/Diff(S 1 ) [27], while if φ is an immersion and φ(S 1 ) is diffeomorphic to S 1 , we obtain the shape space Emb(S 1 , R 2 )/Diff(S 1 ) (roughly, the curve has no self-intersections). There are also smooth (C ∞ ) versions of these, piecewise smooth versions, shapes of bounded variation, and so on. The geometry of these shape spaces is now studied intensively in its own right [24,25] as well as as a setting for computer vision; see [9] for a recent review.
The quotient by Diff(S 1 ) in the shape spaces has the effect of factoring out the dependence on the parameterization. Various methods have been used to achieve this, including: 1. Using a standard parameterization, such as Euclidean arclength. This leaves a dependence on the choice of starting point, i.e., the subgroup of Diff(S 1 ) consisting of translations is not removed.

Moment invariants,
L 0 f (φ(s)) ds, where L is the length of the curve, s is Euclidean arclength and f ranges over a set of basis functions, such as monomials or Fourier modes [1].
3. Currents S 1 φ * α, where α ranges over a set of basis 1-forms on R 2 , such as monomials times dx and monomials times dy [16].
The next step is to consider a Lie group G acting on R 2 which induces an action on the shape space. Here, many methods and types of invariants have been investigated (see e.g. [41]).
The moving frame method is a general approach to constructing invariants. Objects are put into a reference configuration and their resulting coordinates are then invariant. The moving frame or reference configuration method was developing as a way of finding differential invariants (see [31] for an overview). A simple application of the method is the situation considered in Example 1, of n points in the plane under rotations. Consider configurations with z 1 = 0, | arg z 1 | < π. Rotate the configuration so that z 1 lies on the positive real axis. The coordinates of the resulting reference configuration,z 1 z j /|z 1 | 2 , are invariant. In this case, the invariants are well behaved as | arg z 1 | → π, but not as z 1 → 0.
Joint invariants are functions of several points; for example, pairwise distances for a complete set of joint invariants for the action of the Euclidean group on sets of points in the plane.
Differential invariants are functions of the derivatives of a curve at a point. There exist algorithms to generate all differential invariants [31]. For the action of the Euclidean group on planar curves, the Euclidean curvature κ = φ × φ / φ 3 is E(n) invariant (and parameterization invariant), and its derivatives d n κ/ds n with respect to arclength form a complete set of differential invariants.
Semi-differential invariants of a curve are functions of several points and derivatives.
Integral invariants are formed from the moments or the partial moments s s0 f (φ(t)) dt. With some care they can be made parameterization-and basepoint-independent [12,22]. Initially, these invariants appear to have some advantages, being relatively robust and often including some locality. However, they are not always applicable; for example [22], which used regional integral invariants, still requires a point correspondence optimization in order to get a distance between shapes. In addition, groups such as the projective group do not act on any finite subset of the moments [41]. Astrom [6] shows that there are no stable projective invariants for closed planar curves. While local integral invariants are promising, as reported by [17], there may be analytic difficulties in deriving them.
Another example of the moving frame method, which is common in image processing and shape analysis, is centre of mass reduction. The centre of mass of a shape may be moved to the origin in order to remove the translations. This may be calculated by, for example, However, the shapes x = a + sin t, y = 0 have centre of mass equal to 0 for all values of a, even though they are related by translations. Thus, this method would not be robust on any dataset which contained shapes approaching such degenerate shapes. For rotations, the reference configuration method will not be robust if there are objects close to having a discrete rotational symmetry. The underlying problem with the reference configuration approach is that it is attempting to use a set of invariants equal to the dimension of the desired quotient space M/G. This space is almost always non-Euclidean, so such a set can only be found on some subset of M/G of Euclidean topology. The set is then robust only on datasets that are bounded away from the boundary of this subset.
Calabi et al. [10] propose the use of differential invariant signatures for shape analysis, and further argue that these should be approximated in a group-invariant way. For example, for the Euclidean group, the signature is the shape (κ, κ s ) regarded as a subset of R 2 . The claimed advantages of the approach are that the signature determines the shape; that it does not depend on the choice of initial point on the curve or on parameterization by arclength, and that its Ginvariance makes it robust; and that it is based on a general procedure for arbitrary objects and groups. (See [4,10,19,30,34,35,2] for further developments and applications of the invariant signature.) Although the signature at first sight appears to be complete (e.g., Theorem 5.2 in [10]), a more detailed treatment (e.g. [29,18]) highlights the fact that it is not complete on shapes that contain singular parts-straight and circular segments in the Euclidean case. For these parts, the signature reduces to a point. Thus, for shapes that have nearly straight or nearly circular segments, the signature cannot be robust. In addition, while the signature does not depend on the starting point or the parameterization, it takes values in a very complicated set, namely, the planar shapes. To compare the invariants of two shapes requires comparing two shapes. Essentially, the parameterization-dependence has only been deferred to a later stage of the analysis (unless one is content to compare shapes visually).
The claim in [10] that the method's G-invariance makes it robust should be assessed through further analysis and experiment. It would appear to be most relevant for datasets in which the errors due to the presentations of the shapes are comparable to those resulting from the errors in the shapes themselves (i.e., noise) and from the distribution of the objects in a classification problem. Their final point, that the method is extremely general, is a powerful one. Calabi [10] carry out the procedure for the Euclidean and for the 2D affine group. However in §6.7 they report that "the interpolation equations in general are transcendentally nonlinear and do not admit a readily explicit solution," indicating that the method may not succeed for all group actions.
It is also possible to represent the shape as a binary image and apply image-based invariants (which are described next). This has the advantage of working directly in the space R 2 on which the group acts, and avoiding all questions of parameterization, etc., but it does sacrifice a lot of information about the shape. A related approach, which is popular in PDE-based method for curves, is to represent the shape as the level set of a smooth function φ : R 2 → R. This is also parameterization-independent, and retains smoothness, but we have never seen it used in for constructing invariants.

Invariants of images
Most studies of image invariants has been based on grey-level images f : The methods are primarily based on moments or Fourier transforms of the images, and have been highly developed for the translation, Euclidean, and similarity groups, where the linearity of the action and the special structure of these Lie groups means that the approach is particularly fast and robust. The representation theory of the group is of fundamental importance. Abu-Mostafa [1] studies many aspects of moment-based invariants.
We sketch here the Fourier-based invariants for the translation group. Let Ω = [0, 2π] 2 and letf m,n = Ω e −i(mx+ny) f (x, y)dxdy be the Fourier components of the image f . Translations act on each Fourier mode separately, so that the action of translation by (a, b) becomes f m,n → e i(am+bn)f m,n . Early work used the Fourier amplitudes |f m,n | to recognize images up to translation. Clearly, these are ideal from the point of view of computational complexity and of noise, but they are not complete-indeed, the space of N × N -pixel images modulo translations has dimension N 2 − 2, but there are only N 2 /4 independent Fourier amplitudes. The bispectrum [20]f m1,n1fm2,n2f−m1−m2,−n1−n2 , is better, being complete on images in which all Fourier coefficients are non-zero, but it is very large (O(N 4 ))-the problem is similar to that encountered in Example 1. When scalings are included, there are numerical problems related to the truncation of the infinite domain on which the scaling group acts.
Still, the method of Fourier invariants remains a very powerful and effective one, and attempts have been made to extend it to other groups. There is a harmonic analysis for many non-Abelian Lie groups, including, in fact, the Möbius group [36], as well as a general theory for compact non-commutative groups [14]. There are some applications of this theory to image processing [20,40] and to other problems in computational science. [13] discusses a Fourier transform for the hyperbolic group, the 3-dimensional subgroup of the Möbius group that fixes the unit disc. However, the theory appears not to have been developed to the point where it can be used as effectively as the standard Fourier invariants.

Möbius invariants for curves
In this section we consider invariants of such curves for the Möbius group, and derive one suitable for practical computation, demonstrating its application for a set of simple closed curves.

Known Möbius invariants
The most well-known invariant of the Möbius group is the cross-ratio, also known as the wurf, which is based on the ratio of the distances between a set of 4 points: where the invariance means that CR(T z 1 , T z 2 , T z 3 , T z 4 ) = CR(z 1 , z 2 , z 3 , z 4 ) for Möbius transformation T .
Since the Möbius group can send any triple of distinct points in C into any other such triple, there are no joint invariants of 2 or 3 points; the orbits are the configurations consisting of 1, 2, and of 3 distinct points. When n > 4, the set of cross-ratios of all subsets of 4 of the points forms a complete invariant.
For large numbers of points, this set of all cross-ratios has cardinality O(n 4 ) which is impractically large. However, if the dataset of shapes or images is tagged with a small number of clearly-defined landmarks, then some subset of the set of cross-ratios may form a useful invariant. This is the method by which Petukhov [33] was able to identify linear-fractional, Möbius, and projective relationships in biological shapes. For untagged objects, automatic tagging may be possible using critical points (e.g. maxima and minima) of images, and their values; these are homomorphism-and hence Möbius-invariant, and can be identified, even in the presence of noise, by the method of persistent homology [11]. However, such invariants are clearly highly incomplete for shapes and images, and we do not study them further here.
In order to derive differential invariants for the Möbius group, the most useful starting point is the Schwarzian derivative By an abuse of notation, which is standard in the literature (see, for example, the very readable [3]), the same formula (3) is used in three different situations: when z : R → R (used in studying linear fractional mappings in real projective geometry); when z : C → C (used in studying complex analytic mappings); and when z : M → C, M a real 1-dimensional manifold (used in studying the Möbius geometry of curves). We adopt the latter setting so that z is the tangent to the curve. The Schwarzian derivative is then invariant under Möbius transformations φ: and under reparameterizations ψ : M → M transforms as where the last term is the real Schwarzian derivative. Therefore This can be used to construct a distinguished Möbius-invariant parameterization of the curve.
The parameter λ will be chosen so that ImS(z) ≡ 1. This gives The choice ψ (t) = |ImS(z)(t)| achieves this while preserving the sense of the curve, while the choice ψ (t) = − |ImS(z)(t)| achieves it while reversing the sense. Put another way, the is invariant under Möbius transformations and almost invariant under sense-preserving reparameterizations of t, which act as translations in λ (because of the freedom to choose t 0 ). Equivalently, the 1-form dλ = |ImS(z)(t)|dt, known as the Möbius or inversive arclength, is Möbius-and sense-preserving parameterizationinvariant. This is often stated in a form (originally due, according to Ahlfors [3], to Georg Pick) using the Euclidean curvature κ. If the curve is parameterized by Euclidean arclength s, then its tangent θ(s) = z (s)/ z (s) and κ(s) = θ (s). Differentiating again leads to ImS(z)(s) = κ (s), or dλ = |κ (s)|ds.
The 1-form dλ provides a useful discrete invariant, the Möbius length L of the curve: The real part of S(z)(λ) is now a parameterization-invariant Möbius invariant known as the inversive or Möbius curvature [32] where denotes differentiation with respect to arclength s. The set of all differential Möbius shape invariants is then d n κ Möb /dλ n , n ≥ 0. Following Calabi et al. [10], two possible candidate invariants that could be used to recognize Möbius shapes are the function κ Möb (λ) modulo translations and the signature (κ Möb , dκ Möb /dλ)(S 1 ) ⊂ R 2 . These are complete on sections of shapes with no vertices (points where κ (s) = 0). However, since κ Möb requires the 5th derivative of the curve (third derivative of the curvature), it is not robust in the presence of noise, and we do not explore it further.
There are also invariants based on forms of higher degree. We give just one example, the Möbius energy: introduced by O'Hara and Solanes [28]. (See also the related Kerzman-Stein distance [8].) Here du, dv are Euclidean arclength and θ u (resp. θ v ) is the angle between v − u and z (u) (resp. v). It represents a renormalization of the energy of particles on the curve, distributed evenly with respect to Euclidean arclength and interacting under an r −4 potential. (The authors of [28] comment that "due to divergence problems, almost nothing is known about integral geometry [i.e., about invariant differential forms] under the Möbius group".) The Möbius energy has one distinguishing feature compared to the other invariants: its definition depends only on the first derivative of the curve. The singularity at u = v is removable, since the integrand obeys Thus, the energy is defined for C 2 curves, and even near the singularity it depends only on the curvature of the curve.
However, in numerical experiments we found the integrand of (9) difficult to evaluate accurately and invariantly, particularly near the diagonal, while the double integral generated little improvement in robustness. The most negative feature of the energy is its cost: its evaluation apparently requires considering all pairs of points on a discrete curve. There is another reason why invariant 2-forms are less useful than invariant 1-forms like dλ: invariant k-forms distinguish coordinates on M k only up to k-form-(i.e., volume-) preserving maps. These are infinite-dimensional for k > 1, but 1-dimensional for k = 1: the coordinates are determined up to translations. For these reasons, we do not consider the energy, or related invariant 2-forms, any further.
Finally, Möbius transformations map circles to circles, and thus critical points of Euclidean curvature (the 'vertices' of the shape) are Möbius invariants. The four vertex theorem states that all smooth curves have at least four vertices. The number of vertices is a discrete Möbius invariant. Consequently, sections of shapes on which the Euclidean curvature is monotonic, and their Möbius lengths, are Möbius invariant.

Example: Evaluating the Möbius length of an ellipse
Having described a number of possible invariants, and ruled many of them out based on the criteria discussed in Section 2.1, we now provide a concrete example, illustrating and testing some of these constructions on the ellipse z(t) = cos 2πt + 2i sin 2πt, 0 ≤ t ≤ 1. The ellipse is discretized at t i = (i + 1/4)h, i = 0, . . . , n, giving points z i = z(t i ), and the Möbius arclength dλ is calculated in two ways: • From the curvature method, in which the Euclidean curvature κ is calculated in a Euclideaninvariant way 2 by interpolating a circle through 3 adjacent points, and dκ/ds by a finite difference, giving: • From the cross-ratio method, using: Away from vertices (points where dλ = 0), both of these are second-order finite difference approximations to the Möbius length ti+2 ti+1 dλ of the arc z([t i+1 , t i+2 ]), or, dividing by h, to its arclength density dλ/dt. We test their accuracy as a function of the step size h. The total Möbius length of the curve is: L − L h is expected to have dominant contributions of order h 3/2 , due to the singularities at the vertices, and of order h 2 , due to the finite differences used to approximate dλ. In this example, the cross-ratio approximation has errors approximately 0.176 times those of the curvature approximation (see Figure 3). However, due to some pecularities of the singular integrand, the curvature approximation actually gives a slightly more accurate approximation to the Möbius length (see Figure 4). The dominant sources of error can be eliminated by two steps of Richardson extrapolation, first to remove the O(h 3/2 ) error, and then to remove the O(h 2 ) error. This is highly successful and allows the calculation of the Möbius length with an error of less than 10 −10 , even though it is singular and involves a 3rd derivative.
Next, we subject the ellipse to a variety of Möbius transformations. The resulting shapes and errors are shown in Figure 5. The errors increase markedly for the curvature method, which is not Möbius invariant, but are unchanged for the cross-ratio method, which is Möbius invariant. Thus, this experiments supports the argument of [10] that numerical approximations of invariants should themselves be invariant.

The Möbius invariant CR(λ; z, δ)
The method proposed in this paper for closed curves is to parameterize the curve by Möbius arclength, giving z(λ), and to use as an invariant the cross-ratio of all sets of 4 points a distance δ apart. We call this the Shape Cross-Ratio, or SCR.
The shape cross-ratio signature of z is the shape SCR(R; z, δ) ⊂ C.   (10) and (11). The shape cross-ratio is invariant under the Möbius group, and sense-preserving reparameterizations of the curve act as translations in λ. (Reversals of the curve will be considered in §3.5.) The shape cross-ratio signature is invariant under the Möbius group and under reparameterizations of the curve.
For an L-periodic function f : R → C we will denote its Fourier coefficients by: The translation t → t+c acts on the Fourier coefficients as F(f ) n → e −2πicn/L F(f ) n . The Fourier amplitudes |F(f ) n | 2 are invariant under translations, and can be used to recognize functions up to translations, but are clearly not a complete invariant: for a function discretized at N equallyspaced points, and using the DFT, the space of orbits has dimension 2N − 1 and we have only N invariants. The bispectrum is better, being complete on functions all of whose Fourier coefficients are non-zero, but it is a very large set of invariants. Other invariants are F(φ 1 •f ) n F(φ 2 • f ) n for any functions φ 1,2 . Each such choice provides 2N invariants. The choice of φ 1 and φ 2 determines which aspects of f are measured by the invariant. If necessary, several such pairs may be used.

Definition 2
The Fourier cross-ratio of the shape z is: where the Fourier transforms are based on the Möbius length L of z.
In the numerical illustrations we use: and the distance between the invariants of two shapes z 1 and z 2 given by: The motivation here is that the cross-ratio becomes arbitrarily large when two different parts of the curve approach one another. If left untouched (i.e., if we use just φ 1 (w) = w), then these large spikes in SCR(λ; z, δ) will dominate all other contributions to the shape measurement. By scaling them using (13), they will still contribute to the description of the shape, but in a way that is balanced with respect to other parts of the shape. φ 1 and φ 2 take values in the unit disk and φ 2 is sensitive to the main range of features of shapes; only values of SCR near 0 are suppressed, and these are rare.
The invariant SCR(λ; z, δ) is smooth on simple closed curves, and also on most curves with self-intersections (blow-up requires the close approach of two points Möbius distance nδ apart.) It is locally complete, as given z([a, a + 3δ]), the invariant SCR(λ) determines z. We do not know if it is globally complete, i.e., if, given SCR(λ) which is the invariant of some shape, the shape can be determined up to Möbius transformations, because this requires solving a nonlinear functional boundary value problem. Subject to this restriction, the invariant FCR(n; z, δ) is complete except on a residual set of shapes (those for which enough Fourier coefficients in (12) are zero.) We will first study the numerical approximation of SCR and then study its use in recognizing shapes modulo Möbius transformations.
The numerical experiments in Section 3.2 convinced us to approximate the Möbius arclength using the cross-ratio. However, when combined with piecewise linear interpolation to locate points on the curve the required distance δ apart, we found that the resulting values of SCR(λ; z, δ) did not converge as h → 0. This is due to accumulation of errors along the curve, which arise particularly at the vertices due to the singularities there. This prompted us to develop a more refined interpolation method that takes into account the singularities of dλ/dt at the vertices. We call it the modified cross-ratio method: 1. Calculate the square of the Möbius arclength density at the centre of each cell, as If the curve is smooth, this is a smooth function.  |f (τ )| dτ and its inverse, t(λ), used in locating the parameter values at which points a desired length apart are located, exactly (we omit the formulas).
4. The desired points z(t(λ + nδ)) are calculated using linear interpolation from the known values z(ih).

5.
The cross-ratio is evaluated at N points equally spaced in λ, giving SCR(iL/N ; z, δ) for i = 1, . . . , N , and the Fourier invariant FCR evaluated using two FFTs.
The resulting cross-ratio is globally second-order accurate in h. Its accuracy could be increased for smooth curves using higher order interpolation, but the calculation of the inverse t(λ) would be much more complicated and the method would be less robust. The errors in the length L of the ellipse used in Section 3.2 as calculated by this method is shown in Figure 4. It behaves extremely reliably over a wide range of scales of h, and its error after one Richardson extrapolation is observed to be O(h 4 ), which indicates that the singularities at the vertices have been completely removed.
The paramater δ is the length scale on which SCR(λ; z, δ) describes the shape. However, if δ → 0 then SCR(λ; z, δ) → κ Möb ±i: the real part becomes extremely nonrobust and the imaginary part yields no information [32]. In the experiments in this paper we have used δ = L/8. Choosing L/δ ∈ Z seems to yield somewhat improved accuracy as the same values of z are used repeatedly.
As the method relies on parameterization by Möbius arclength, it is unavoidably sensitive to noise in the data. We test its sensitivity for a noise model in which each point on the discrete curve is subject to normally distributed noise of standard deviation ε. The dependence of the error on and N is shown for length and cross-ratio invariants in Figure 6. Clearly both are sensitive to relatively small amounts of noise. However, some positive features can also be seen: (i) The errors in L and FCR(n; z, δ) are both O(h 2 ) in the absence of noise.
(ii) The errors in FCR(n; z, δ) are much smaller than those in L-their relative errors are about 8 times smaller. (For the ellipse, L ≈ 6.86 and SCR(·; z, δ) 2 ≈ 0.65.) (iii) The error in FCR(n; z, δ) appears to saturate at about 13% as ε increases and as N increases.
Point (iii) is particularly striking. It appears to hold because (a) δ is chosen to be proportional to L, and thus when the noise is high, the chosen points stay roughly in their correct places; and (b) noise in the chosen points is averaged out by the Fourier transform, which remains dominated by its first few terms. This effect is illustrated in Figure 7 in which a single noise realization is illustrated for value of ε from 0 to 10 −2 . Even though L is overestimated by a factor of 100, the signature cross-ratio SCR(λ; z, δ) is still recognizable.

Comparison with shape registration
We now compare the results of the Möbius invariant FCR(n; z, δ) with direct registration of shapes. Given two shapes z and w we define the G-registration of z onto w as: Different choices of norm in (15) will give different registrations; we have used the H 1 norm where the constant α was chosen as 0.1, a value which made both contributions to the norm roughly equal. One of the peculiarities of the Möbius group is that z may register very well onto w while w registers poorly onto z. This happens when z has a distinguished feature which can be squashed, thus minimizing its contribution to r G (z, w). Therefore in our experiments we use the 'distance': Note that this should not be regarded as any kind of 'ground truth' for the G-similarity of z and w. It is not G-or parameterization-invariant. However, as we shall see it does correspond remarkably well to the Möbius invariants described earlier.
To calculate d G (z, w) numerically, we discretize Diff + (S 1 ) by piecewise linear increasing functions with 16 control points and perform the optimization using Matlab's lsqnonlin, with initial guesses for ϕ chosen to be each of 4 rotations and 3 scale factors for ϕ, and initial ψ chosen to be the identity.   Figure 9: Scatter plot of distance between all 120 pairs of 16 shapes with respect to (i) the H 1 distance between the shapes, as defined in (15)- (17), and (ii) the distance between their Möbius invariants, defined in (12)- (14). The correlation coefficient is 0.85. The registration between the 4 pairs marked with circles is illustrated in Figures 10-13.  . This pair has the 2nd closest Möbius invariants (distance 0.0697), see (14)), is 3rd closest after Möbius registration, and 42rd closest after similarity registration.  . This pair has the 3rd closest Möbius invariants (distance 0.0798), is 9th closest after Möbius registration, and 96th closest after similarity registration.
The shapes are shown in Figure 8. They are all simple closed curves which we take to be positively oriented. The 120 pairs of distinct shapes are registered in both directions, and the scatter plot between this distance and the 2-norm of the distance between their invariants is shown in Figure 9. The correlation (0.85) is extremely striking, and suggests that this pair of measures may be related by bounded distortion (1), implying that they meet many of the requirements that we identified in Section 2.1; they are also quick to compute and provide a good numerical approximation.
Some examples of the registrations in the similarity and Möbius groups are shown for four close pairs in Figures 10-13. Including the invariants L and V in the list of invariants did not improve the correlation. Note that the errors in FCR(·; z, δ) observed in Figure 4 are small enough to allow the separation of all but the closest pairs of shapes, regardless of the level of noise ε.

Reversals and reflections
Orientation reversals and reflections are examples of actions of discrete groups and can, in theory, be handled by any of the approaches in Section 2.2. First, consider the sense-reversing reparameterizations, z(t) → z(−t). These map the SCR invariant as SCR(λ; z, δ) → SCR(−λ; z, δ) and hence FCR(n) → FCR(−n). It is convenient to pass to the equivalent invariants: for which x n is invariant for n ≤ 0 and x n → −x −n for n > 0. Suppose that we wish to identity curves with their reversals, that is, to work with unoriented shapes. Some options are the following: 4. As the group action is of a standard type, one can work in unreduced coordinates (x n ) together with a natural metric induced by the quotient, such as some function of x − y and the Fubini-Study metric on projective space, cos −1 (|x T y|/ x y ).
5. Finally, and most easily in this case, for finite groups one can represent points in the quotient as entire group orbits. A suitable metric is then where G is the group. Although this is impractical for large finite groups, here G is Z 2 .
Similar considerations apply to recognizing shapes modulo the full inversive group, generated by the Möbius group together with a reflection, say z →z. The reflection maps FCR(n) → FCR(n) and hence is an action of the same type as reversal (a sign change in some components). The reflection symmetry of the ellipse in Figure 2, for example, can be detected by the reflection symmetry of the cross-ratio signature SCR(λ) in Figure 7. (Its second discrete symmetry, a rotation by π, is manifested in Figure 7 by the signature curve retracting itself twice.) The invariants developed here are for closed shapes. They can be adapted for other types of shapes (for example, shapes with the topology of two disjoint circles), but as the topology gets more complicated (for example, shapes consisting of many curve segments) the problem becomes significantly more difficult.

Möbius invariants of images
Let f : C → [0, 1] be a smooth grey-scale image. Diffeomorphisms act on images by ϕ·f := f •ϕ −1 . It is easy, in principle, to adapt the Möbius shape invariant SCR to images by computing level sets of f , each of which is an invariant shape for which SCR can be calculated. In addition, if ϕ is conformal, the orthogonal trajectories of the level sets, i.e., the shapes tangent to ∇f , are also invariant shapes. In the neighbourhood of a simple closed level set, coordinates (λ, µ) can be introduced, where λ is Möbius arclength along the level set and µ is Möbius arclength along the orthogonal trajectories. The quantity CR(z(λ, µ), z(λ, µ + δ), z(λ + δ, µ + δ), z(λ + δ, µ)), calculated from the cross-ratio of 4 points in a square is then invariant under the Möbius group, and reparameterizations of the level set act as translations in λ.
In practice, however, the domain of this invariant is quite restricted. The topology of level sets is typically very complicated and the domain of f may be restricted, so that level sets can stop at the edge of the image. Restricting to level sets of grey-scales near the maximum and minimum of f helps, but this is a severe restriction. Instead, we shall show that the extra information provided by an image, as opposed to that provided by a shape, determines a differential invariant signature using only 3rd derivatives, compared to the 5th derivatives needed for differential invariants of shapes. Because of this, we do not develop the cross-ratio invariant (18) any further here.
Proposition 1 Let f : C → R be a smooth grey-scale image. Let R ⊂ C be the regular points of f . On R, define n := ∇f ∇f , is a subset of R 3 that is invariant under the action of the Möbius group on images.
Proof 1 As defined above, n is the unit vector field normal to the level sets of f . Let n ⊥ be the unit vector tangent to the level sets given by n ⊥ i = ε ij n j . From the Frenet-Serret relation n s = κn ⊥ , where s is arclength along the level sets, we have that the curvature of the level sets is κ = n ⊥ · n s = n ⊥ · ((n ⊥ · ∇)n) = ε ij n j ε kl n l n i,k = (δ ik δ jl − δ il δ jk )n j n l n i,k = n j n j n k,k − n k n i n i,k = n k,k (because n j n j = 1 ⇒ n i n i,k = 0 for all k) = ∇ · n (20) Recall that the Möbius arclength of the level sets of f is dλ := |κ s |ds. Under any conformal map, the scaling along and normal to the level sets is the same, and thus ds and 1/ ∇f both scale by the same factor. Therefore |κ s |/ ∇f is invariant, as is its square |κ s |/ ∇f 2 . The sign of κ s is also invariant under Möbius transformations, resulting in the given invariant λ t = κ s / ∇f 2 .
The invariant λ n arises in the same way from the orthogonal trajectories, whose curvature is ∇ · n ⊥ = ∇ × n.
Example 2 As a test image we take the smooth function and calculate its invariant signature before and after the Möbius transformation with parameters a = 0.9 + 0.1i, b = 0.1, c = 0.1 + 0.4i, d = 1.
on the domain [−1, 1] 2 . The invariants are approximated by finite differences with mesh spacing 1/80, corresponding to 161 × 161 pixel images. The resulting signature surface, shown in figure  17 in R 3 is quite complicated, but in figures 14 and 15 we show contour plots for a the function and the Möbius transformed version in the (f, λ n ) and (f, λ t ) coordinate planes respectively. Figure 17 shows the invariant signature of this same function f , and it can be seen that the signatures are similar.   show just the signature curve corresponding to the level set f −1 (0.5). This depends only on the first 3 derivatives of f on the level set. Because λ n and λ t take values in [−∞, ∞], we use coordinates (arctan(λ t /4), arctan(λ n /4)). Clearly, even this very limited portion of the signature serves to distinguish the Möbius-related pairs extremely sensitively. In some cases, the invariants change extremely rapidly along the level set, so that even though they are evaluated accurately, the resulting contours of the Möbius-related pairs do not overlap. This would need to be taken into account in the development of a distance measure on the invariant signatures.

Example 4
In this example we illustrate the extreme sensitivity of the invariant signature by evaluating it on 9 very similar images, together with their Möbius transformations. Each original image is a blob function generated as in Example 3, but with parameters varying only by ±5%. The Möbius transformations have the form 1/(1 + cz) where c is normally distributed with standard deviation 0.1. The 0.5-level contours of the original and transformed images are shown in Figure 20, and their signatures in Figure 21. The signature is extremely sensitive to tiny changes in the image, but not to Möbius transformations.
We do not have a full understanding of the properties of this invariant signature with respect to the criteria listed in Section 2. It is certainly fast, small, local, and lacks redundancy and suppression. It has a good numerical approximation on smooth (or smoothed) images. It is complete? That is, given an image, does its signature surface determine the image up to Möbius? Suppose we are given a small piece of signature surface, parameterized by (u, v), say. We are given are given three functionsf (u, v),λ n (u, v), andλ t (u, v), and need to determine (by solving three PDEs) three functions f (x, y) (the image), x(u, v), and y(u, v) (the coordinates). Typically the solution of these PDEs will be determined by some boundary data. This suggests that distinct   Figure 18 left (shown in blue) and for the corresponding images in Figure 18 right (shown in red). The domains are [−π/2, π/2] 2 . The signature curves distinguish the Möbius-related pairs very sensitively; only tiny finite difference errors are visible. However, some errors related to insufficient resolution of the signature curves are clearly visible. images with the same signature are parameterized by functions of 1 variable; a kind of near completeness that may be good enough in practice. Although very sensitive, the fact that it is not continuous at critical points means that it does not have good discrimination in the sense of Section 2. (It falls into the 'more false negatives' region of Fig. 1.) Near nondegenerate critical points, the signature blows up in a well-defined way, so it is possible that there exists a metric on signatures which leads to robustness and good discrimination.

Conclusion
In this paper we have developed Möbius invariants of both curves and images, and proposed computational methods to evaluate both, demonstrating them on a variety of examples. In Section 2 we identified a set of properties that are important for invariants, principally that there was a small set of invariants that were quick to compute, numerically stable, robust (so that noisy versions of the same curve have similar invariants) and yet sufficiently discriminatory (so that different objects have different invariants).
While differential invariants are not generally robust when dealing with noise, they offer good discrimination and are cheap to compute; this lead us to the Möbius arc-length. The cross-ratio is more robust, but requires a large set of points to be evaluated, and blows up as the pairs of points approach each other. In order to make this reparameterization-invariant, we used a Fourier transform. This lead to a method of computing Möbius invariants that satisfies the properties that we have outlined, as is demonstrated in the numerical experiments, see particularly Figure 9.
For images, the extra information means that it is possible to compute a relatively simple three dimensional signature based on the function value at each point together with two functions of the Möbius arclength, along and perpendicular to level sets of the image intensity. It is computationally cheap, extremely sensitive to non-Möbius changes in the image, but insensitive to Möbius transformations of the image.