Currents and finite elements as tools for shape space

TL;DR

This paper introduces a new approach for representing and analyzing shapes using currents and finite element methods, enabling robust numerical discretization and shape comparison in nonlinear shape spaces.

Contribution

It develops the theory of currents for shape spaces, analyzes their properties, and provides a finite element discretization for practical shape analysis.

Findings

Current map determines shapes under certain conditions.

Finite element discretization is effective for shape computation.

Approach is robust to noise and suitable for numerical implementation.

Abstract

The nonlinear spaces of shapes (unparameterized immersed curves or submanifolds) are of interest for many applications in image analysis, such as the identification of shapes that are similar modulo the action of some group. In this paper we study a general representation of shapes that is based on linear spaces and is suitable for numerical discretization, being robust to noise. We develop the theory of currents for shape spaces by considering both the analytic and numerical aspects of the problem. In particular, we study the analytical properties of the current map and the $H^{-s}$ norm that it induces on shapes. We determine the conditions under which the current determines the shape. We then provide a finite element discretization of the currents that is a practical computational tool for shapes. Finally, we demonstrate this approach on a variety of examples.

Equations

$$[\varphi ](\alpha ):={\int }_M{\varphi }^{\ast }\alpha$$

$$[\varphi ]:{\Lambda }^m(N)\to \mathbb{R},[\varphi ](\alpha ):={\int }_M{\varphi }^{\ast }\alpha ={\int }_{\varphi (M)}\alpha .$$

$$[\cdot ]:\mathrm{Imm}(M,N)\to {\Lambda }^m(N{)}^{\ast }.$$

$${\int }_M(\varphi \circ \psi {)}^{\ast }\alpha ={\int }_M{\psi }^{\ast }{\varphi }^{\ast }\alpha ={\int }_{\psi (M)}{\varphi }^{\ast }\alpha ={\int }_M{\varphi }^{\ast }\alpha$$

$$(\alpha ,\beta {)}_{H_{\sigma }^s}:={\int }_{\Omega }\sum _{\genfrac{}{}{0ex}{}{i=1,2}{0\le |k|\le s}}{\sigma }^{2|k|}(\genfrac{}{}{0ex}{}{s}{k})(D^k{\alpha }_i)(D^k{\beta }_i)\text{d}{x}_1\text{d}{x}_2$$

$$\varphi =\left({\varphi }_1,{\varphi }_2\right),$$

$$\left[\varphi {|}_V\right](\alpha )={\int }_V({\alpha }_1\circ \varphi \cdot {\varphi }_1^{\prime }+{\alpha }_2\circ \varphi \cdot {\varphi }_2^{\prime })\text{d}t.$$

$$\left|\left[\varphi {|}_V\right](\alpha )\right|\le (b-a)\left({\Vert {\alpha }_1\Vert }_{L^2(\varphi (S^1))}{\Vert {\varphi }_1\Vert }_{W^{1,\infty }(V)}+{\Vert {\alpha }_2\Vert }_{L^2(\varphi (S^1))}{\Vert {\varphi }_2\Vert }_{W^{1,\infty }(V)}\right)=C<\infty .$$

$$\Vert [\varphi ]{\Vert }_{H^{-s}{\Lambda }^1(\Omega )}:=\underset{\genfrac{}{}{0ex}{}{\alpha \in H^s{\Lambda }^1(\Omega )}{\Vert \alpha {\Vert }_{H^s}=1}}{\sup }[\varphi ](\alpha )$$

$$d({\varphi }_1,{\varphi }_2):=\Vert [{\varphi }_1]-[{\varphi }_2]{\Vert }_{-s}:=\Vert [{\varphi }_1]-[{\varphi }_2]{\Vert }_{H^{-s}{\Lambda }^1(\Omega )}.$$

$$A:=\left|[\varphi ](\alpha )-[\psi ](\alpha )\right|.$$

$$A=\left|{\int }_{\varphi (S^1)}\alpha -{\int }_{\psi (S^1)}\alpha \right|.$$

$$A={\int }_R|\text{d}\alpha |.$$

$${\left({\int }_{R}fg\text{d}x\wedge \text{d}y\right)}^2\le \left({\int }_R{f}^2\text{d}x\wedge \text{d}y\right)\left({\int }_R{g}^2\text{d}x\wedge \text{d}y\right)$$

$$A$$

$$\le \Vert \alpha {\Vert }_{H^1}{\left(area(R)\right)}^{\frac{1}{2}}$$

$$\le \Vert \alpha {\Vert }_{H^1}(C\delta {)}^{\frac{1}{2}}$$

$$\Vert [\varphi ]-[\psi ]{\Vert }_{H^{-1}}=\underset{\Vert \alpha {\Vert }_{H^1}=1}{\sup }\left|[\varphi ](\alpha )-[\psi ](\alpha )\right|\le (C\delta {)}^{\frac{1}{2}},$$

$$A\le \sum _{i=0}^n{d}_i|f|\text{d}x\wedge \text{d}y$$

$$A\le \Vert \alpha {\Vert }_{H^1}{\left(\sum _{i=0}^n{d}_i^{2}area(R_i)\right)}^{\frac{1}{2}}$$

$${\Delta }_{\epsilon }$$

$$={\int }_{\tilde{\xi }(\epsilon ,S^1)}{\pi }^{\ast }\alpha -{\int }_{\tilde{\xi }(0,S^1)}{\pi }^{\ast }\alpha$$

$$={\int }_{{\Gamma }_{\epsilon }}{\pi }^{\ast }\alpha$$

$$={\int }_{M_{\epsilon }}\mathrm{d}{\pi }^{\ast }\alpha .$$

$$\displaystyle\frac{d}{d\varepsilon}[\xi(\varepsilon,\cdot)](\alpha)\big{|}_{\varepsilon=0}$$

$$={\int }_{\tilde{\xi }(0,S^1)}{\mathrm{i}}_{{\tilde{\xi }}_{\ast }\frac{\partial }{\partial u}}\mathrm{d}{\pi }^{\ast }\alpha$$

$$={\int }_{\varphi (S^1)}{\mathrm{i}}_X\mathrm{d}\alpha .$$

$$[\varphi ](\alpha )=(\beta ,\alpha {)}_{H^s}$$

$$\Vert [\varphi {]}_{-s}{\Vert }^2:=\Vert \beta {\Vert }_{H^s}^2=[\varphi ](\beta ).$$

$$(\alpha ,\beta {)}_{H_{\sigma }^1}$$

$$\displaystyle=\int_{\Omega}\Big{(}\sum_{i=1}^{2}\alpha_{i}(1-\sigma^{2}\nabla^{2})\beta_{i}\Big{)}\,\mathrm{d}x_{1}\mathrm{d}x_{2}+\int_{\partial\Omega}\beta_{in}\alpha_{i}\mathrm{d}s$$

$$\langle\delta_{\phi(S^{1})},f\rangle=\int_{0}^{1}f\big{(}\phi(t)\big{)}\|\phi^{\prime}(t)\|\,\mathrm{d}t.$$

$$(\alpha )$$

$$={\int }_0^1⟨\alpha (\varphi (t)),{\varphi }^{\prime }(t)⟩\mathrm{d}t$$

$$={\int }_0^1⟨\alpha (\varphi (t)),\mathbf{t}(t)⟩\Vert {\varphi }^{\prime }(t)\Vert \mathrm{d}t$$

Full text

Keywords: Currents, finite elements, shape space, image analysis.

MSC 2010: 32U40m 62M40, 65D18, 74S05

1 Introduction

“Shape,” wrote David Mumford in his 2002 ICM address Mumford (2002), “is the ultimate nonlinear thing.” The set of smooth simple closed curves is an example of a shape space; in the study of shape one seeks analytic and numerical methods to work with such curves: to compare, to classify, to recognise, to evolve—to understand them. As is usual in mathematics, nonlinear things are constructed out of linear things, namely vector spaces, using simple operations: maps, quotients, and open subsets.

The aforementioned curves are constructed by first considering parameterised curves $\phi\colon S^{1}\to{\mathbb{R}}^{2}$, and then identifying curves that differ only by a reparameterisation. The set of smooth parameterised curves consists of all smooth embeddings $\mathrm{Emb}(S^{1},{\mathbb{R}}^{2})$ of $S^{1}$ into ${\mathbb{R}}^{2}$. A smooth reparameterisation corresponds to composition from the right by a diffeomorphism of $S^{1}$, i.e., by an element of the set $\mathrm{Diff}(S^{1})$. Thus, the set of (images of) simple closed curves is given by the quotient $\mathrm{Emb}(S^{1},{\mathbb{R}}^{2})/\mathrm{Diff}(S^{1})$. More generally, let $M$ and $N$ be manifolds of dimension $m$ and $n$ respectively. A shape space is the set of $m$-dimensional submanifolds of $N$ that are diffeomorphic to $M$, which is realized as $\mathrm{Emb}(M,N)/\mathrm{Diff}(M)$. (This is a nonlinear analogue of the Grassmannian of linear subspaces of a specific dimension.) There are other shape spaces as well, such as the immersed shapes $\mathrm{Imm}(M,N)/\mathrm{Diff}(M)$, shapes in which $M$ has a boundary, piecewise smooth shapes, and oriented shapes.

It is desirable to recognise examples of the same shape, those that are identical up to the action of $\mathrm{Diff}(M)$ and possibly noise or other obfuscation. This can be studied by finding a metric that is blind to changes inside the class, or by finding a representation of the shapes that is invariant under the action of $\mathrm{Diff}(M)$. An example of such a representation is the differential invariant signature of shapes, an influential new paradigm introduced by Calabi et al. (1998). For example, to recognize planar curves up to Euclidean transformations a signature curve is $(\kappa,\kappa_{s})(S^{1})\subset{\mathbb{R}}^{2}$, where $\kappa$ is Euclidean curvature and $\kappa_{s}$ its derivative with respect to arclength. This is clearly invariant under Euclidean motions and under reparameterisations. As Calabi et al. comment, “The recognition problem includes a comparison principle that would be able to tell whether two signature curves are close in some sense. Thus, we effectively reduce the group-invariant recognition problem to the problem of imposing a ‘metric’ on the space of shapes but now by ‘shape’ we mean the signature curve, not the original object.” Such metrics can be constructed through shape currents.

The method of shape currents, first suggested by Glaunès et al. (2008), is based on embedding the nonlinear shape space in a vector space endowed with a metric, thereby allowing the construction of flexible families of metrics on shape spaces that are easy to compute. It is robust to noise and provides for control of the resolution and accuracy of the representation.

In previous work, metrics on shape currents have been combined with optimization routines for registration, typically by deforming the shape by left action of a diffeomorphism group (see Section 1.1 below for more details of the use of currents in shape analysis). Here, we take the viewpoint of Calabi et al. that left symmetries (registration) are taken care of by computing a signature curve, so the only remaining step is to impose a metric on shape space. This way we obtain a distance on shapes modulo any classical transformation group; see Example 10.

The central object in the paper is the current map (see 1), which takes a function $\phi\colon M\to N$ and associates it with

[TABLE]

where $\alpha\in\Lambda^{m}(N)$ is an $m$-form on $N$. The current map is invariant under orientation-preserving reparameterizations (see 1). Continuing from Glaunès et al. (2008), our aim is, in broad terms, to study this map and the shape distances it induces. Let us summarize the main results.

We fix a class of functions $\phi$ and a topology on them, namely that of the Lipschitz immersions $\phi\colon S^{1}\to\Omega\subset{\mathbb{R}}^{2}$, and likewise fix a class of forms $\alpha$ and a topology on them. Since $[\phi]$ is a linear map from 1-forms to the reals, it is natural to first adopt a Sobolev norm on 1-forms and then to demand that $[\phi]$ be a continuous map, so that we can adopt the corresponding operator norm on currents.

In Section 2, we show the following results (see 2 and 3 for the definitions of $H^{s}\Lambda^{1}(\Omega)$, $H^{-s}\Lambda^{1}(\Omega)$ and $\mathrm{LipImm}(S^{1},\Omega)$):

• •

The map $\alpha\mapsto[\phi](\alpha)$ defined in (1) is continuous and linear in $\alpha$ (2), thus showing that $[\phi]$ lies in $H^{-s}\Lambda^{1}(\Omega)$. This generalizes (Glaunès et al., 2008, Proposition 1) to Lipschitz curves.

• •

The current map from $\mathrm{LipImm}(S^{1},\Omega)$ to $H^{-s}\Lambda^{1}(\Omega)$ is Hölder-continuous for $s\geq 1$ (3)

• •

The current map (with same domain and codomain as above) is differentiable for $s\geq 2$ (4). This regularity is important: to use invariants to recognize shapes, they must have the property that nearby shapes have nearby invariants.

• •

For a given subspace of the space of functions $S^{1}\to\Omega$, and a given linear space of 1-forms on $\Omega$, we show that the current map $\phi\mapsto[\phi]$ essentially determines the shape of $\phi$ (5).

In the numerical part (Section 3), we discretize our construction. The currents can be evaluated for all $\alpha$ in a finite element space, which gives a flexible and general representation of shapes. Specifically, we introduce spaces of finite elements $V$ on $M$ and $W$ on $\Lambda^{m}(N)$ and evaluate the current $[\phi_{V}]|_{W}$, where $\phi_{V}$ is the approximation of $\phi$ in $V$. A simple example is shown in Figure 1. In this case, piecewise constant elements determine a piecewise constant approximant that interpolates the shape at the element edges. The norm on currents restricts to $W$ in a natural way, yielding a discretized norm on shapes.

In Section 3, we show the following results:

• •

In Section 3.2, we demonstrate that the quadrature errors in evaluating the currents are typically small and that the method is robust in the presence of noise. This robustness stems from the cancellation property of oscillatory integrals. The currents can be accurately computed even for very rough shapes (not even Lipschitz) and for noisy shapes.

• •

In Section 3.3, we study how accurately the discretized currents determine the shapes. This is a question in approximation theory.

• •

In Propositions 8 and 9 we show that the order of approximation is 2 for piecewise constant, 3 for piecewise linear, and 5 for piecewise quadratic elements in $W$.

• •

In Section 3.4, we introduce a discretization of the metric on shapes, so that each shape is approximated by a point in a vector space equipped with an Euclidean metric.

Finally, we give several numerical examples of the geometry of shape space induced by the discretization in Section 4. The method does not only compare pairs of shapes, it provides a direct approximation of shape space and its geometry: we present numerical experiments (e.g., Examples 9 and 13) applying Principal Components Analysis directly to the current representation in order to successfully separate classes of shapes.

1.1 Related Literature

Currents have already been used in shape analysis, primarily for curve and surface matching. In Vaillant and Glaunès (2005) a surface in 3D was represented with currents defined on a surface mesh, and a norm was computed to enable the matching of two surfaces using the currents. A framework is developed to allow shape currents to be matched in the spirit of the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework Beg et al. (2005), and the method was demonstrated on surfaces representing shapes and hippocampi. A variation on this approach for curves rather than surfaces was developed in Glaunès et al. (2008) (where the currents are referred to as vector-valued measures). Again, the aim is a matching algorithm where curves can be deformed onto each other.

In Durrleman et al. (2009) a different benefit of the linear representation provided by currents is recognised, which is that they provide a useful space in which to perform statistical analysis of the deformations between shapes. This was originally considered in Glaunès and Joshi (2006), but there is a difficulty that the mean (template) shape has to be defined in such a way that both the shape and deformation are amenable to statistical analysis. A Matching Pursuit algorithm is defined in Durrleman et al. (2009) to provide a computationally tractable representation of shape currents as the number of shapes grows.

The use of currents for matching was extended by Charon and Trouvé in Charon and Trouvé (2014) as functional currents, where a function is added to the current representation so that the deformation of a shape and some function defined on its surface can be considered simultaneously. The same authors also considered how to deal with cases where the orientation property of currents is undesirable. The fact that a large spike in the appearance of the shape cancels out in the current representation as the positive and negative contributions are virtually identical is a benefit when considering the currents for dealing with noisy representations of shapes. However, in cases where these spikes can truly exist, or where there is orientation information in the image, but it can differ arbitrarily by sign, such as in Diffusion Tensor Images of the brain, the oriented manifold is a disadvantage. This leads to the consideration of varifolds in Charon and Trouvé (2013), where the registration of some directed surfaces is demonstrated.

While there are many other computational approaches to shape space, as far as we are aware they all involve determining a point correspondence between shapes via optimization and/or working directly in the nonlinear shape space; see for example Celledoni et al. (2016); Bauer et al. (2014).

In contrast to the aforementioned work, we focus here solely on shape currents as a way to induce distances and compute statistics on $\mathrm{Emb}(M,N)/\mathrm{Diff}(M)$; we assume that registration (left matching) has already been taken care of, for example through the signature curve, as suggested by Calabi et al. (1998).

2 Currents and their induced metric on shapes

For any natural number $m$, we denote by $\Lambda^{m}(N)$ the space of smooth $m$-forms on a manifold $N$. An $m$-current is an element of $\Lambda^{m}(N)^{*}$, the topological dual of $\Lambda^{m}(N)$. Later, in Section 2.2, we will introduce another topology on $\Lambda^{m}(N)$, and we will still call a current an element of the corresponding dual.

Currents were introduced by de Rham De Rham (1973). They are natural generalisations (or completions) of the pairing by integration between an $m$ dimensional submanifold $M\subset N$ and smooth $m$–forms on $N$. Currents are instrumental in geometric measure theory, where they are used to study a very wide class of (not necessarily smooth) subsets of ${\mathbb{R}}^{n}$, for example in minimal surface problems Morgan (2009). In that field, the functions $\phi$ are typically Lipschitz and the differential forms are smooth. We allow nonsmooth shapes such as those represented by Lipschitz functions, but as we want to discretize the differential forms by finite elements, we will allow the forms to be nonsmooth as well.

2.1 De Rham currents and their invariance

Our parameterization invariant map (1) is motivated by the signed area $\int_{\phi(S^{1})}y\,\textrm{d}x$ enclosed by a closed curve. This is induced from a 1-form ($y\,\textrm{d}x$) on ${\mathbb{R}}^{2}$.111This also happen to be Euclidean-invariant, but this is not relevant to the sequel. Currents are invariant under orientation-preserving diffeomorphisms of $M$—sense-preserving reparameterizations in the case of curves—and hence can be used to factor out the group $\mathrm{Diff}^{+}(M)$, the subgroup of orientation-preserving diffeomorphisms. As can be seen by considering the length and area of a noisy curve (see the example in Figure 7), currents are very robust to noise.

Although we will mostly focus in the sequel on the case of $M=S^{1}$ and $N=\Omega$, where $\Omega$ is a domain in ${\mathbb{R}}^{2}$, we formulate the definition of the current map in a more general setting:

Definition 1

Let $M$ and $N$ be oriented manifolds of dimensions $m$ and $n$, with $m\leq n$. We denote by $\mathrm{Imm}(M,N)$ the space of immersions from $M$ to $N$. Let $\phi\colon M\to N$ be an immersion. We define $[\phi]$ to be the linear function on forms given by

[TABLE]

The current map is the corresponding map from parameterized manifolds to currents:

[TABLE]

Proposition 1

The current map defined in 1 is invariant with respect to orientation-preserving reparameterizations, that is $[\phi\circ\psi]=[\phi]$ for all orientation-preserving diffeomorphisms $\psi$ of $M$. In other words, the current map $\phi\mapsto[\phi]$ is $\mathrm{Diff}^{+}(M)$-invariant (where $\mathrm{Diff}^{+}(M)$ denotes the space of orientation-preserving diffeomorphisms on $M$), and it induces a map from $\mathrm{Imm}(M,N)/\mathrm{Diff}^{+}(M)$ to $\Lambda^{m}(N)^{*}$.

Proof

This is just the statement that integration of forms is well-defined, i.e., independent of the choice of coordinates. Specifically,

[TABLE]

for orientation-preserving $\psi$. ∎

Currents measure some particular aspects of curves. Consider the case of planar curves, and in particular the curve $\phi(t)=(\cos(t),0)$ for $0\leq t\leq 2\pi$. This shape retraces itself in opposite directions, so $[\phi](\alpha)=0$ for all $\alpha$; currents cannot distinguish this curve from the 0 curve.

In the case that $\phi(M)$ is a submanifold of $N$, $[\phi]$ is the $m$-current of integration on $\phi(M)$.

2.2 Properties of the current map

Before choosing a particular norm in the space of forms, let us motivate our choice. For the current map to be defined (1), we need to be able to take traces of a form on any curves (see 1). In particular, we want this operation to be continuous (see 2). This is related to the definition of a reproducing kernel Hilbert space, in which traces on points are supposed to exist. Notice however that our setting is strictly more general, as, for instance, $H^{1}({\mathbb{R}}^{2})$ is not a reproducing kernel Hilbert space, but traces on curves are well defined.

So, for its metric and computational aspects, we adopt the following $H^{s}$ Sobolev inner product on 1-forms:

Definition 2

Given a non-negative integer $s$, a scale parameter $\sigma$, and a domain $\Omega$ of ${\mathbb{R}}^{n}$, we denote by $H_{\sigma}^{s}\Lambda^{1}(\Omega)$ the Hilbert space of forms with scalar product defined by

[TABLE]

where $\alpha$ and $\beta$ are 1-forms with coordinates $\alpha=\sum_{i}\alpha_{i}\textrm{d}x_{i}$ and $\beta=\sum_{i}\beta_{i}\textrm{d}x_{i}$. The dual of that space will be called $H_{\sigma}^{-s}$ currents, and denoted by $H_{\sigma}^{-s}\Lambda^{1}(\Omega)$. We will often omit the scale parameter $\sigma$ in the sequel for the sake of readability.

To understand why $\sigma$ is called a scale parameter, suppose that the domain $\Omega={\mathbb{R}}^{n}$ for the moment. Then to any positive scalar $\lambda$ we can associate the scaling $\lambda\cdot x:=\lambda x$. This function $(\lambda\cdot)$ acts on 1-forms by pull-back, that is (with a slight abuse of notation) $\lambda\cdot\alpha:=(\lambda^{-1}\cdot)^{*}\alpha$. The resulting form is $\lambda\cdot\alpha(x)=\lambda^{-1}\sum_{i=1,2}\alpha_{i}(\lambda^{-1}x_{1},\lambda^{-1}x_{2})\textrm{d}x_{i}$, so we see that $(\lambda\cdot\alpha,\lambda\cdot\beta)_{H^{s}_{\sigma}}=\lambda^{-1}(\alpha,\beta)_{H^{s}_{\sigma/\lambda}}$. We will also see in Section 2.3 that $\sigma$ determines the length scale at which distances between shapes are be measured.

The properties of the map from curves to $H^{-s}$ currents depends not only on the topology of the currents, but also on that of the curves. We wish to allow very large classes of curves—indeed, one of the strengths of currents is that they do allow this. We thus define the following space of curves.

Definition 3

We define the space of Lipschitz immersions, that is, curves $\phi\colon S^{1}\to\Omega$ such that the components of $\phi$ are Lipschitz and the tangent vector $\phi^{\prime}$, wherever it is defined (which is almost everywhere), is nonzero. We call the space of such curves $\mathrm{LipImm}(S^{1},\Omega)$.

Remark 1

Given a bounded domain $\Omega$ in ${\mathbb{R}}^{n}$, with Lipschitz boundary $\partial\Omega$, a typical function $f$ in $H^{1}(\Omega)$ is not continuous and is only defined almost everywhere in $\Omega$. Moreover, $\partial\Omega$ has $n$-dimensional Lebesgue measure zero; hence there is no direct meaning we can give to the expression “$f$ restricted to $\partial\Omega$”. The notion of a trace operator and the trace theorem resolves this issue for us Adams and Fournier (2003). More generally, passing from functions in $H^{s}(\Omega)$ to their traces on surfaces of codimension 1 results in a loss of smoothness corresponding to half a derivative.

The following result generalizes (Glaunès et al., 2008, Proposition 1), from piecewise $C^{1}$ to Lipschitz immersions.

Proposition 2

Let $s\geq 1$ and let $\phi\colon S^{1}\to\Omega$ be a Lipschitz immersion. Then $[\phi]\colon H^{s}\Lambda^{1}(\Omega)\to{\mathbb{R}}$ is a bounded linear operator.

Proof

Each Lipschitz immersion $\phi$ can be written as:

[TABLE]

where $\phi_{1}$ and $\phi_{2}$ are Lipschitz functions from $S^{1}$ to $\mathbb{R}$.

Parameterizing $S^{1}$ by $t\in[0,1)=I$ and letting $\left\{U_{i}\right\}$ denote a finite system of open sets covering $I$, functions $\phi_{1}$ and $\phi_{2}$ (5) are Lipschitz functions from $U_{i}$ to $\mathbb{R}$. We recall Evans (1998) that $\phi_{1}|_{U}$ and $\phi_{2}|_{U}$ are Lipschitz if and only if they belong to $W^{1,\infty}(U)$, the space of functions with essentially bounded first weak derivative.

Let $(a,b)=V\in\left\{U_{i}\right\}$ and $\alpha=\alpha_{1}\,\textrm{d}x+\alpha_{2}\,\textrm{d}y\in H^{s}\Lambda^{1}({\mathbb{R}}^{2})$. In coordinates we have:

[TABLE]

Since $\phi_{1}|_{V}$ and $\phi_{2}|_{V}$ both belong to $W^{1,\infty}(V)$, and $\alpha_{1}\circ\phi$, $\alpha_{2}\circ\phi\in L^{2}(\phi(V))$, by the trace theorem (see 1), we obtain:

[TABLE]

Since $I$ is compact we can find a constant $\tilde{C}$ which works for all sets $U_{i}$, so $[\phi]$ is bounded. ∎

2 implies that the dual (operator) norm

[TABLE]

is well defined. That is, $[\phi]$ is an element of the Sobolev dual $H^{-s}\Lambda^{1}(\Omega)$. We measure the similarity of shapes by their distance in this dual: for two shapes $\phi_{1}$, $\phi_{2}$, their distance is (see Figure 2 for a pictorial version of this):

[TABLE]

To sum up, currents map shapes into a Hilbert space (in a highly nonlinear way) and we measure the distance between shapes using the norm on that Hilbert space. At first sight this appears to be a wasteful representation, as it uses functions on the higher dimensional space $\Omega\subset{\mathbb{R}}^{2}$ instead of on $S^{1}$. However, in finite dimensional examples in which quotient spaces are represented using invariants, it is common that large numbers of invariants are needed222For example, when $S^{1}$ acts on ${\mathbb{C}}^{n}$ by $z_{i}\mapsto e^{i\theta z_{i}}$, the set of invariants $\bar{z}_{i}z_{j}$, $1\leq i,j\leq n$—$n^{2}$ real invariants in all—is complete, and $n^{2}$ is much larger than $\mathrm{dim}({\mathbb{C}}^{n}/S^{1})=2n-1$. One can find smaller complete sets, limited by the dimension of the smallest Euclidean space into which ${\mathbb{C}}^{n}/S^{1}$ can be embedded. However, in such sets the individual invariants are more complicated Bandeira et al. (2014), so that the total complexity is $\mathcal{O}(n^{2})$.. Furthermore, this approach allows one to represent much larger classes of objects than the smooth embeddings, such as weighted, nonsmooth, and immersed shapes.

As remarked earlier, it is vital that the map from curves to $H^{-s}$ currents be continuous, so that nearby shapes have nearby $H^{-s}$ currents. The topology of Lipschitz immersions is the direct product of $W^{1,\infty}$ in each component. That is, two Lipschitz immersions $\phi$, $\psi$ are close if the suprema of $|\phi_{1}(t)-\psi_{1}(t)|$, $|\phi_{2}(t)-\psi_{2}(t)|$, $|\phi^{\prime}_{1}(t)-\psi^{\prime}_{1}(t)|$, and $|\phi^{\prime}_{2}(t)-\psi^{\prime}_{2}(t)|$ over $0\leq t<1$ are all small.

Proposition 3

Let $s\geq 1$. Then the current map $[\cdot]\colon\mathrm{LipImm}(S^{1},\Omega)\to H^{-s}\Lambda^{1}(\Omega)$ is Hölder-continuous with exponent $1/2$; in particular, it is continuous.

Proof

It suffices to prove the proposition for $s=1$, for if two linear functions are close in $H^{-1}$ then they are close in $H^{-s}$ for $s>1$. First consider the case of embeddings. Let $\phi$ be a fixed Lipschitz embedding and let $\psi$ be any Lipschitz embedding with $\|\phi_{1}-\psi_{1}\|_{W^{1,\infty}(S^{1})}<\delta$, $\|\phi_{2}-\psi_{2}\|_{W^{1,\infty}(S^{1})}<\delta$. Let $\alpha\in H^{1}\Lambda^{1}(\Omega)$. We need to estimate

[TABLE]

First, we have:

[TABLE]

Let $R$ be the region enclosed between $\phi(S^{1})$ and $\psi(S^{1})$ (see Figure 3). From Stokes’s theorem:

[TABLE]

In coordinates, $\alpha=\alpha_{1}\textrm{d}x+\alpha_{2}\textrm{d}y$, which gives $\textrm{d}\alpha=(\partial_{x}\alpha_{2}-\partial_{y}\alpha_{1})\textrm{d}x\wedge\textrm{d}y$. Defining $f=\partial_{x}\alpha_{2}-\partial_{y}\alpha_{1}$, we have $\textrm{d}\alpha=f\textrm{d}x\wedge\textrm{d}y$. Now, applying the Cauchy-Schwarz inequality

[TABLE]

with $g\equiv 1$, we get

[TABLE]

where $C$ is a constant depending on $\phi$ (approximately equal to the length of $\phi$). Therefore, for all such $\psi$ we have

[TABLE]

establishing the claim.

For immersions that are not embeddings, Eq. (6) is modified to take into account any intersections. Let the region $R$ between $\phi$ and $\psi$ be $R=\cup_{i=0}^{n}R_{i}$ where $R_{0}$ is the nonoverlapping part and $R_{1},\dots,R_{n}$ are the overlapping parts. Then

[TABLE]

where $d_{0}=1$ and each $d_{i}$ is either 0 or 2, depending on the orientation of the boundary curves of each $R_{i}$ (see Figure 3). Applying the Cauchy–Schwarz inequality with $g=\sum_{i=0}^{n}d_{i}\mathcal{X}(R_{i})$ gives

[TABLE]

As the number of intersections is fixed by the choice of $\phi$, again we have $A\leq(C\delta)^{\frac{1}{2}}$, establishing the claim. ∎

We now show that if the forms are smooth enough ($s\geq 2$), then the current map is differentiable. In the case $s=1$, 7 indicates the map $[\cdot]$ is then not differentiable.

Proposition 4

Let $s\geq 2$. Then the current map $[\cdot]\colon\mathrm{LipImm}(S^{1},\Omega)\to H^{-s}\Lambda^{1}(\Omega)$ is differentiable.

Proof

Let $\xi(u,\cdot)$ be a $C^{1}$ curve in $\mathrm{LipImm}(S^{1},\Omega)$ with $\xi(0)=\phi$ and $\partial_{u}|_{u=0}\xi(u)=X$. Let $\alpha\in H^{s}\Lambda^{1}(\Omega)$. Recall that for a 2-form $\omega$, the interior derivative $\mathrm{i}_{X}\omega$ is defined by $\mathrm{i}_{X}(Y):=\omega(X,Y)$. We will show that the derivative of the current map at $\phi$ in direction $X$ is $\int_{\phi(S^{1})}{\rm i}_{X}{\rm d}\alpha$. By the trace theorem (see 1), ${\rm d}\alpha\circ\phi\in H^{s-2}(S^{1})$, and hence because $s\geq 2$, the derivative exists.

To establish this formula for the derivative, we first consider the case that the curves $\xi(u,\cdot)$ are embeddings. Let $\varepsilon>0$. Let $\tilde{\xi}(u,t):=(u,\xi(u,t))$ and $\pi(u,\tilde{\xi})=\tilde{\xi}$. Then $\tilde{\xi}([0,\varepsilon],S^{1})$ is a tubular surface $M_{\varepsilon}$ with boundary $\Gamma_{\varepsilon}:=\tilde{\xi}(\varepsilon,S^{1})-\varepsilon(0,S^{1})$. Then, using Stokes’s theorem in the last line,

[TABLE]

Letting $\varepsilon\to 0$,

[TABLE]

In the case that the curves $\xi(u,\cdot)$ are immersions, the range of the curves may be lifted from ${\mathbb{R}}^{2}$ to ${\mathbb{R}}^{3}$ and the curves perturbed slightly at the crossings so that they become embeddings in ${\mathbb{R}}^{3}$. Then the same formula for the derivative holds in ${\mathbb{R}}^{3}$ and, letting the perturbation tend to zero, in ${\mathbb{R}}^{2}$.

∎

Thus when $s\geq 2$ one can define a continuous Riemannian metric on shapes as the restriction of $H^{-s}\Lambda^{1}(\Omega)$ to the currents of shapes. While many families of Riemannian metrics on shapes have been studied Michor and Mumford (2007), this one appears to be new. However, in this paper we do not use the induced Riemannian metric, but rather the (‘straight line’) subset metric illustrated in Figure 2, which is far easier to compute.

2.3 Representers of shapes

From 2 the Riesz representation theorem applies: the current $[\phi]$ determines a unique $\beta\in H^{s}\Lambda^{m}(\Omega)$ called the (Riesz) representer of $[\phi]$, which satisfies:

[TABLE]

for all $\alpha\in H^{s}\Lambda^{m}(\Omega)$. We can also write

[TABLE]

That is, two shapes are close in $H^{-s}$ if the representers of their currents are close in $H^{s}$.

From the definition of the representer, which is a PDE in weak form, the representer satisfies an elliptic PDE with source concentrated on the shape. We now take a closer look at the PDE in strong form in the case where both the representer and the shape are smooth.

Suppose now that $\phi(S^{1})$ is smooth, and that the representer $\beta$ is smooth. From Green’s theorem, for $s=1$ we have

[TABLE]

where $\mathrm{d}s$ is arclength of the boundary, and $\beta_{in}$ is the derivative of $\beta_{i}$ normal to the boundary of $\phi(S^{1})$.

Now let $\delta_{\phi(S^{1})}$ be the measure on $\Omega$ with support $\phi(S^{1})$ that obeys

[TABLE]

for all smooth $f\colon\Omega\to{\mathbb{R}}$. Using this definition, we have

[TABLE]

where $\langle\alpha,\mathbf{t}\rangle$ denotes (with a slight abuse of notation) any extension of that function outside the curve $\phi(S^{1})$. Therefore

[TABLE]

for all $\alpha\in H^{1}\Lambda^{1}(\Omega)$. This gives the Helmholtz equation

[TABLE]

For a small scale $\sigma$, the representer $\beta$ can be identified with a vector field concentrated near the curve and roughly tangent to it (an example is shown in Figure 10). The representer can be written in terms of the Green’s function of the Helmholtz operator. For example, when $\Omega={\mathbb{R}}^{2}$, the Green’s function of $(1-\sigma^{2}\nabla^{2})$ is

[TABLE]

where $\mathsf{K}_{0}$ is the modified Bessel function of the 2nd kind; then

[TABLE]

As remarked earlier, for immersions, the currents do not determine the shape, as parts of the shapes that retrace themselves are invisible to currents. For embeddings the situation is better:

Proposition 5

Let $\phi_{1},\phi_{2}\colon S^{1}\to\Omega$ be two Lipschitz embeddings of $S^{1}$ and let $s=1$. If $[\phi_{1}]=[\phi_{2}]$ then the curves represent the same oriented shape.

Proof

Suppose that a point $x$ in $\phi_{1}(S^{1})$ is not in $\phi_{2}(S^{1})$. Then there is a neighbourhood of $x$ which does not intersect $\phi_{2}(S^{1})$ either. Choose a 1-form $\alpha$ that has support in this neighbourhood. As $\phi_{1}(S^{1})$ is an embedding, the form $\alpha$ can be chosen such that $[\phi_{1}](\alpha)\neq 0$. However, $[\phi_{2}](\alpha)=0$, which gives a contradiction. ∎

However, currents do determine the shape of immersions if there is some control over the self-intersections:

Proposition 6

Let $\phi_{1},\phi_{2}\colon S^{1}\to\Omega$ be two Lipschitz immersions, each with a finite number of self-intersections such that at each self-intersection the tangent vectors are continuous and distinct. Let $s=1$. If $[\phi_{1}]=[\phi_{2}]$ then the curves represent the same oriented shape.

Proof

As in the proof of 5, the currents determine the images of $\phi_{1}$ and $\phi_{2}$ away from self-intersections. At the self-intersections, the hypothesis on continuity and distinctness of the tangent vectors allows the non-self-intersecting pieces of the shapes to be joined together in a unique way, thus determining the same oriented shape. ∎

The assumption on the self-intersections is necessary. Even without parts of curves that retrace themselves, self-intersections—either with equal or discontinuous tangent vectors—prevent the current from recognising the shape up to $\mathrm{Diff}(S^{1})$ reparameterizations (see Fig. 4). The current sees only the image of the curve.

Insight into the $H^{-s}$ shape metric is obtained by considering the target domain ${\mathbb{R}}^{2}$ and shapes consisting of vertical lines with periodic boundary conditions in $y$. 7 illustrates both the non-differentiability of the current map $[\cdot]$ in the $s=1$ case and the difference between the $s=1$ and $s=2$ metrics. The $s=1$ metric weights nearby portions of the shapes more heavily than the $s=2$ metric does. It also shows the role of the length scale $\sigma$; roughly, all curves more than Euclidean distance $\sigma$ apart are an equal distance apart in the $H^{-s}$ metrics.

We first prove an elementary Lemma.

Lemma 1

Suppose that $H$ is a reproducing kernel Hilbert space of functions on ${\mathbb{R}}$. Let $\delta_{x}\in H^{*}$ denote the evaluation at the point $x\in{\mathbb{R}}$. Suppose further that the kernel is translation invariant, i.e., the representer for $\delta_{x}$ takes the form $x^{\prime}\mapsto K(x^{\prime}-x)$. Then the distance between $\delta_{0}$ and $\delta_{\epsilon}$ is

[TABLE]

Proof

In general, if $K_{x}$ denotes the representer of $\delta_{x}$, we have

[TABLE]

Using the translation invariance of the kernel, we have $K_{\epsilon}(\epsilon)=K_{0}(0)$, which finishes the proof. ∎

Assuming periodic boundary conditions in $y$, the representer of a vertical line at $x=x_{0}$ is $K((x-x_{0})/\sigma)$ where $K$ is the corresponding one-dimensional kernel. We define the distance per unit length to be the distance in one dimension, with the corresponding one-dimensional kernel.

Proposition 7

The distance per unit length between two straight lines a (Euclidean) distance $\epsilon$ apart is $\mathcal{O}(\epsilon^{1/2})$ when $s=1$ and $\mathcal{O}(\epsilon)$ when $s=2$.

Proof

Using 1, the distance by unit length is thus $\sqrt{2(K(0)-K(\epsilon))}$.

For $s=1$, $K(x)=\frac{1}{2}\mathrm{e}^{-|x|}$ and the distance per unit length between the lines is $(1-\mathrm{e}^{-\epsilon/\sigma})^{1/2}=\mathcal{O}(\epsilon^{1/2})$ as $\epsilon\to 0$.

For $s=2$, $K(x)=\frac{1}{4}\mathrm{e}^{-|x|}(1+|x|)$ and the distance per unit length between the lines is $(\frac{1}{2}(1-\mathrm{e}^{-(\epsilon/\sigma)}(1+\epsilon/\sigma)))^{1/2}=\mathcal{O}(\epsilon)$ as $\epsilon\to 0$. ∎

We now turn to a practical implementation of shape representation using currents that is based on the finite element method.

3 Discretization of currents by finite elements

The discretization of the representation of shapes by currents proceeds in two steps:

• (i)

Discretization of the space of currents and approximation of the currents of shapes; and

• (ii)

Discretization of the metric on currents.

We shall consider these separately, as their numerical properties are somewhat independent; item (i) determines how accurately the shapes themselves are represented, while item (ii) determines the geometry of the induced discrete shape space.

3.1 Discretization of currents

We first return to the general setting in which we work with immersions $\phi\colon M\to N$ where $M$ and $N$ are manifolds and with their currents $[\phi](\alpha)=\int_{M}\phi^{*}\alpha$ for $\alpha\in H^{s}\Lambda^{m}(N)$. The discretization of currents requires the choice of three things:

• (i)

A space $V$ of finite elements on $M$;

• (ii)

A space $W$ of finite elements on $\Lambda^{m}(N)$;

• (iii)

A method of evaluating or approximating $[\phi_{V}]|_{W}$, where $\phi_{V}$ is the finite element representation of $\phi$ in $V$.

We will explore analytically and numerically the ability of $[\phi_{V}]|_{W}$ to represent shapes in different settings.

To illustrate the ideas, we first consider an extremely simple example, namely sets of $n$ points on a line.

Example 1

Let $M=\{1,\dots,n\}$ and $N={\mathbb{R}}$. A shape is then an unordered set of $n$ points in ${\mathbb{R}}$ with isotropy the group of permutations. In this case $M$ does not need to be discretized. We choose a uniform mesh with spacing $\Delta x$ on ${\mathbb{R}}$ and let $W$ be the piecewise polynomials of degree at most $d$ on the mesh. Let $w(x)\mathrm{d}x\in W$. Then

[TABLE]

The simplest case is $d=0$, i.e., piecewise constant elements. For these, the currents count how many points are in each cell. Therefore these elements represent the shape with an accuracy of $\mathcal{O}(\Delta x)$, as there is no way to tell where in each cell the points are located.

The next case is $d=1$, i.e., piecewise linear elements. The piecewise constants determine how many points are in each cell, and the linear elements determine $\sum_{j\colon x_{j}\in W_{i}}x_{j}$, that is, the mean location of the points in cell $W_{i}$. If there is at most 1 point in each cell, then the representation is perfect, as the points are located exactly. If there is more than 1 point in a cell, then these elements represent the shape with an accuracy of $\mathcal{O}(\Delta x)$. Notice how the use of currents factors out the parameterization of the shape, i.e., it is invariant under permutations of $M$.

With piecewise elements of degree $d$, the number of points in each cell and their first $d$ moments are determined; if there are at most $d$ points in each cell, the points are located exactly. This is because the moment equations on $[0,\Delta x]$ (for example) are $\sum_{j\colon x_{j}\in W_{i}}(x_{j})^{i}=c_{i}$. For given values of the moments $c_{i}$, this set of polynomial equations has total degree $d!$, hence at most $d!$ solutions; but if there is one real solution, then any permutation of the $x_{j}$ is a solution; hence the moments determine the points up to ordering.

We now consider our main example of oriented closed planar curves. Let $M=S^{1}$ and let $N=\Omega$, a domain in ${\mathbb{R}}^{2}$. We will choose $V$ to be the continuous piecewise polynomials of a given degree on a uniform mesh on $S^{1}$ and $W$ to be either the discontinuous or the continuous piecewise polynomials of a given degree on a fixed mesh on $\Omega$. We will take the finite element representation $\phi_{V}$ of $\phi$ to be the element of $V$ that interpolates $\phi$ at the finite element nodes. The current $[\phi_{V}]|_{W}$ is then given by the integral of a piecewise polynomial function. This can be evaluated exactly; however, in this study we evaluate $[\phi_{V}]_{W}$ by quadrature, either the midpoint rule for piecewise linear elements or Simpson’s rule for piecewise quadratics. As we are interested in fairly compact representations of shapes, we will pick the meshsize of $V$ to be much smaller than the meshsize of $W$.

In the following subsections we study the behaviour of finite element currents with regard to (i) quadrature errors, with and without noise; (ii) accuracy of representation of shapes; and (iii) accuracy of the induced metric on shapes.

3.2 Quadrature errors and robustness of currents.

One strong motivation for considering currents, as opposed to other possible shape invariants such as those based on arclength (cf. 1), is that—essentially because they are signed—currents are expected to be robust against noise and to function well on quite rough shapes. We present three examples measuring different aspects of robustness.

Example 2

In this example the shapes are rough, but there is no noise.

Figure 5 shows the quadrature errors for 3 rough shapes as a function of the number of points used to discretize them. In this and the following example, the $H^{-1}$ shape norm is discretized using polynomials of degree 9 on $[-1,1]^{2}$, although only the quadrature error is reported. Shapes (a) and (b) are continuous, but not Lipschitz continuous, and yet the quadrature errors are well under control, and for shape (b) are even $\mathcal{O}((\Delta s)^{2})$ (where $\Delta s$ is the mean mesh spacing).

Example 3

In this example, shown in Figure 6, the shape is smooth, but different levels of noise are added. Thus, the exact value to which we compare is the zero noise, zero mesh spacing limit. Two different quadratures (the midpoint and Simpson’s rule) are compared. Although high levels of noise can dominate the quadrature error, they do not prevent the computation of highly accurate values for the currents until the noise level is actually greater than the mesh spacing $\Delta s$.

Example 4

We compare the sensitivity of currents to that of arclength-based currents in Figure 7. The reference shape is a unit line segment, discretized with uniform mesh spacing $\Delta s$, and independent normally-distributed noise of mean 0 and standard deviation $\epsilon$ is added to the $x$ and $y$ components of each point except the endpoints. The current $\int y\,{\textrm{d}}x$ is computed by the midpoint approximation $\sum(y_{i+1}-y_{i})(x_{i+1}-x_{i})$. For $\epsilon\ll\Delta s$, the errors accumulate like a sum of random variables, the error being $N(0,\epsilon(\Delta s)^{\frac{1}{2}})$. For $\epsilon\gtrsim\Delta s$, nonlinear effects create a larger error of size $\mathcal{O}(\epsilon(\Delta s)^{-\frac{1}{2}})$. Overall, good results are obtained even with noise levels $\epsilon=\mathcal{O}(\Delta s)$. In contrast, the positive nature of length means that the arclength approximation $\sum((x_{i+1}-x_{i})^{2}+(y_{i+1}-y_{i})^{2})^{\frac{1}{2}}$ accumulates errors of size $\mathcal{O}((1+(\epsilon/\Delta s)^{2})^{\frac{1}{2}}-1)$. These are $\mathcal{O}(\epsilon^{2}/(\Delta s)^{2})$ for small $\epsilon$ and $\mathcal{O}(\epsilon/\Delta s)$ for larger $\epsilon$. Overall, $\epsilon=o(\Delta s)$ is necessary for acceptable results.

3.3 Accuracy of representation of shapes

In this section we consider simple closed planar shapes and study how accurately they are represented by discontinuous finite elements. Let $M=S^{1}$ and let $N=\Omega$, a domain in ${\mathbb{R}}^{2}$. Let $\phi\colon S^{1}\to\Omega$ be a smooth embedding. Take a triangular mesh on $\Omega$ consisting of triangles $\mathcal{T}$ of maximum diameter $\Delta x$. Let $W$ be the space of discontinuous piecewise polynomial finite elements of degree $\leq d$. We are studying the effect of the discretization of the ambient space $N$, so we do not discretize $M$; we assume that all currents are evaluated exactly. We assume that the currents $[\phi]|_{W}$ are given and we want to know how accurately they determine the shape.

First we consider piecewise constant elements, that is, $d=0$.

Proposition 8

For smooth embeddings $\phi\colon S^{1}\to\Omega\subset{\mathbb{R}}^{2}$, and a triangular mesh $\mathcal{T}$ of sufficiently small mesh size $\Delta x$, the currents of $\phi$ evaluated on 1-forms constant on each triangle $T\in{\mathcal{T}}$ determine a piecewise linear approximation $\hat{\phi}$ of $\phi$ of pointwise second order accuracy.

Proof

The currents are

[TABLE]

If the mesh is sufficiently fine, then these currents are either zero (if the curve does not intersect $T$), or they record the jumps in $x$ and $y$ of that part of the curve that lies in $T$. The elements on which the currents are nonzero therefore determine the set of elements whose interiors intersect the curve $\phi$. If the mesh is sufficiently fine and the curve is in general position, then these elements form a ‘discrete topological circle’, a set of triangles each sharing an edge with exactly two others (see Figure 9).

We now consider finding a continuous shape $\hat{\phi}$ with the same currents as $\phi$. The currents have two degrees of freedom per triangle, as do shapes that are linear on each triangle; we therefor seek such an approximant $\hat{\phi}$. For any values of the currents of $\phi$ on piecewise constants, i.e., for any values of the jumps in $x$ and $y$, there are at most two line segments in $T$ with endpoints on the edges of $T$ whose currents take on these values. (See Figure 8.) If there are two such line segments, then they join different pairs of edges. The line segment that joins two edges that are part of the known discrete topological circle can then be chosen. See Figure 9.

This piecewise linear approximation $\hat{\phi}$ to $\phi$ interpolates $\phi$ at the edges of the elements, and, as $\phi$ is assumed to be smooth, obeys

[TABLE]

on each triangle $T$. It therefore determines the shape to second order accuracy. ∎

Next we consider the improvement that can be obtained using discontinuous piecewise linear or quadratic elements, i.e., $d=1$, 2.

Proposition 9

Let $\mathcal{T}$ be a triangular planar mesh. Let

[TABLE]

Then for sufficiently smooth $\phi\colon S^{1}\to{\mathbb{R}}^{2}$ in general position and for $\mathcal{T}$ sufficiently fine, the integrals of the first $k$ 1-forms in $V_{x}$ and $V_{y}$, $2\leq k\leq 5$, over $\phi(S^{1})\cap T$ for each $T\in\mathcal{T}$ determine an $\mathcal{O}((\Delta x)^{k})$-accurate approximation of $\phi(S^{1})$.

Proof

From Proposition 8, the piecewise constant currents determine the occupied triangles. Suppose that the shape can be written on an occupied triangle $T$ in the form $y=g(x)$. If the derivatives of $g$ are bounded we use the 1-forms in $V_{x}$ only. Otherwise, the curve can be written in the form $x=f(y)$ where $f$ has bounded derivatives and we use the 1-forms in $V_{y}$ only. Without loss of generality, we consider the first case.

The case $k=2$ is covered in Proposition 8.

For $k=3$, we know the intersection points of the curve with the edges of each occupied triangle and, in addition, the value of the current $\int_{U}y\,\textrm{d}x$ for each $U=\phi(S^{1})\cap T$, $T\in\mathcal{T}$. We take coordinates in which the shape is $y=g(x)$ and the $x$-intersections are at $x=0$ and $x=h$, where $h=\mathcal{O}(\Delta x)$; this fixes the parameterization. We consider the quadratic approximation $y=\hat{g}(x)$ to $g(x)$ that interpolates $g(x)$ at $x=0$ and at $x=h$ and has the same value of $\int_{0}^{h}y\,\textrm{d}x$. This yields the approximation

[TABLE]

where

[TABLE]

At $x=sh$, $0\leq s\leq 1$, the approximation error is

[TABLE]

That is, piecewise linear elements determine the shape to third order accuracy.

For $k=4$, we know, in addition, the current $\int_{U}xy\,\textrm{d}x$, and we choose a cubic approximation. There is a unique such cubic, and it yields a 4th order approximation to $g(sh)$ with leading order error

[TABLE]

For $k=5$, we know, in addition, the current $\int_{U}y^{2}\,\textrm{d}x$, and we seek a degree 4 approximant. The equations for the coefficients of the approximant are now nonlinear, of total degree 2. They have two real solutions. One has leading order error less than $0.00003h^{5}|g^{(5)}(0)|$, and the other has leading order error less than $0.003h^{5}|g^{(5)}(0)|$. This establishes 5th order accuracy for $k=5$. ∎

Note that in practice, the integrals of both $V_{x}$ and $V_{y}$ would be used, but on most triangles they do not provide independent information.

We briefly discuss several things we learn from this proposition:

The currents determine smooth shapes very accurately on sufficiently fine meshes. The errors (less than $0.008h^{3}|g^{(3)}|$ for $k=3$, $0.0006h^{4}|g^{(4)}|$ for $k=4$, and $0.00003h^{5}|g^{(5)}|$ for $k=5$), are less than twice that of Chebyshev interpolation. 2. 2.

If the finite element space $W$ has $k$ degrees of freedom per triangle (here $k=\mathcal{O}(d^{2})$), the order of approximation appears to be $\mathcal{O}(k)$. That is, the higher dimensionality of the target manifold $N$ does not seem to be important. 3. 3.

The inherent nonlinearity of the approximation need not be an obstacle. Despite the nonlinearity of the approximation, leading to quadratic equations when $k=5$, the best approximant can be chosen systematically. Nevertheless, we anticipate that at very high degrees $d$ the nonlinearity may render it difficult to reconstruct an accurate approximation. 4. 4.

The situation here is an example of a moment problem. If we choose coordinates in which the base of the triangle $T$ is at $y=0$, and consider the domain:

[TABLE]

then we are being given the moments $\iint_{E_{i}}x^{m}y^{n}\,\textrm{d}x\textrm{d}y$ (for some set of values of $m,n$) and are asking how accurately we can reconstruct $g(x)$. When $n>1$ we have a nonlinear approximation problem about which, as far as we know, little is known. 5. 5.

The moment problem considered in this section for discontinuous piecewise polynomial finite elements is nonlinear, but at least it is entirely local. We anticipate that analogous results to Proposition 9 hold for the approximation of shapes by the currents of continuous piecewise polynomial finite elements.

3.4 Discretization of the metric on shapes

Recall that the norm of the shape metric is defined by Eqs. (7,8):

[TABLE]

If $W\subset H^{s}\Lambda^{1}(\Omega)$ then both $[\phi]$ and the $H^{s}$ inner product may be restricted to $W$ to yield a finite element approximation of the representer $\beta$ and a metric on $W^{*}$. That is:

[TABLE]

In coordinates, let $w_{1},\dots,w_{k}$ be a basis of $W$ and let ${\mathbf{b}}$, ${\mathbf{f}}$ be the coordinate vectors of $\beta$ and $[\phi]$, respectively; that is, $f_{i}=[\phi](w_{i})$ and $\beta=\sum_{i=1}^{k}b_{i}w_{i}$. Then

[TABLE]

where

[TABLE]

is the matrix of the $H^{s}$-metric restricted to $W$. For $s=1$, $G$ is a linear combination of the mass and stiffness matrices of $W$. Note that (13), (14) amount to a standard finite element solution of the inhomogeneous Helmholtz equation (11). Then we have

[TABLE]

In practice, we use the Cholesky decomposition of $G$ to represent $\beta$ in an orthonormal basis, i.e., we compute $\tilde{\mathbf{b}}=G^{\frac{1}{2}}{\mathbf{b}}$ so that

[TABLE]

In this way each shape maps to a point in a standard Euclidean space ${\mathbb{R}}^{|W|}$ and standard techniques such as Principal Components Analysis can be applied. We think of $\|\cdot\|_{W^{*}}$ as providing a highly compressed or approximate geometric representation of shape space.

Although the choice $s=0$ does not make sense at the continuous level—the ‘representer’ would be a delta function supported on the shape, which is not in $L^{2}$—it does make sense at the discrete level. An example is provided by the piecewise constant finite elements considered in Section 3.3. The representer is nonzero on the triangles that intersect the shape; it is similar to a discrete greyscale drawing of the shape. As the currents of piecewise constant finite elements know the intersection of the shape with the edges of the mesh, they already provide a very sensitive discretization of the metric on shapes. However, it does not converge as $\Delta x\to 0$. Moreover, in this approximate metric, all shapes that intersect with different edges are seen as equally far away. For example, the computed distance between two straight segments located at $x=0$ and $x=d$ is proportional to

[TABLE]

For these reasons we do not consider discontinuous elements further.

An example of an $H^{1}$ representer is shown in Figure 10, calculated using a single square element together with polynomial currents of degree less than 10. The Helmholtz equation (11) smears out the shape over a length scale $\sigma$ ($1/\sqrt{10}$ in this example). This allows shapes to be sensitive to each other’s positions over lengths of order $\sigma$, which is typically many times the mesh spacing.

For any $s$, one possible choice of $W$ is to take polynomials up to some degree $N$. This gives a spectral method. However, we do not expect spectral accuracy because the representers are not smooth, they are only in $H^{s}$. As in finite element solutions of PDEs, the mesh size and order $p$ of the elements can be adjusted depending on the application. But unlike PDEs, in most shape applications it is not necessary to have small errors with respect to the continuous problem. Rather, each choice of $W$ provides a different approximation or description of shape space. In some applications $W$ may be relatively low-dimensional.

For $s=1$, we can take any triangulation of $\Omega$ with $W$ the continuous piecewise polynomials of degree $d\geq 1$; these lie in $H^{1}$ so the above construction applies directly. For $s>1$ we use the same elements, and form the same metric $G$ associated with $s=1$, but determine ${\mathbf{b}}$ by solving

[TABLE]

Such representers satisfy higher-order Helmholtz equations with slightly different boundary conditions than those implied by Eqs. (7) and (8). In practice, this difference is immaterial as long as the curve is located away from the boundaries, since the representer tends to zero far from the curve. Moreover, this approximation error is outweighed by the ease of solving (15) and of using standard finite elements.

4 Examples

We have implemented this approach to currents using finite elements in Python using FEniCS and Dolfin Logg et al. (2012); Alnaes et al. (2015). The implementation is available at Benn et al. (2018). We define a rectangular mesh with continuous Galerkin elements of predefined order. To compute the invariants, we use the tree.compute_entity_collisions() function to identify which mesh cells intersect with the curve. We then evaluate the basis functions of those cells and update the computation of the currents as the sum of the basis elements of intersecting mesh elements weighted by the size of the cell. Solving for the representer in the relevant norm is based on a matrix solve $\mathbf{G}b_{x}=f_{x}$, where $\mathbf{G}$ is the tensor representation of the norm (see equation (14)), $b_{x}$ is a function on the finite element space, and $f_{x}$ is the vector of $x$ currents (i.e., an integral of $w_{i}(x)\,\textrm{d}x$ over the curve), and similarly for $b_{y}$ and $f_{y}$. We used a length scale of $\sigma=1/\sqrt{10}$ throughout, so the $H^{-1}$ norm that we used was computed as m = 1./10inner(grad(u),grad(v))dx () + uvdx (), where u and v were trial and test functions on the function space, respectively. In order to compute the $H^{-2}$ norm it is necessary to perform a second solve, based on the matrix built from xvdx ().

We present a series of experiments demonstrating the use of our implementation on a variety of test cases. The first few examples are chosen to demonstrate the robustness of the approach. They demonstrate the reparameterization of the shape, its numerical stability of the method in both norms, and robustness with curves that lose differentiability, for example at corners. These examples are followed by a set of experiments demonstrating that the representation of the curves that are produced is consistent with human perceptions of the shapes.

Example 5

In the first experiment we next demonstrate that the representation in the $H^{-1}$ and $H^{-2}$ norms is invariant to reparameterization of the curves. The curve is a bowtie shape discretized with 512 points. The equally-spaced positions of these points were then perturbed by a Gaussian random variable of standard deviation 0.1 (recall that the radius of the circle was 0.5) and the points resorted into monotonically increasing order of arclength. Figure 11 shows the representers for the unperturbed shapes on the left, and the perturbed ones on the right using piecewise linear finite elements and a meshsize of $10\times 10$. Visually, there seems to be no difference between them. The difference in the computed norms between the original shape and the reparameterized one was of the order of $10^{-4}$ for both metrics.

Example 6

We examine the convergence of the discretized metric with respect to the meshsize and the order of the finite elements. A circle of radius 0.5 in the domain $\Omega=[-1,1]^{2}$ discretized with 5000 equally spaced points, so that quadrature errors are negligible. The grid is a uniform triangular mesh on an $M\times M$ grid. Using elements of order 1 (piecewise linear) up to order 4, we computed the difference in each of the $H^{1}$ and $H^{2}$ norms as the meshsize $M$ was increased in powers of 2 from 1 up to 128, i.e., we compute $\|[\phi]\|_{H^{-s},2M}-\|[\phi]\|_{H^{-s},M}$. The results are shown in Figure 12. Reference lines illustrating convergence of order $1$ for $H^{-1}$ and order 2.5 for $H^{-2}$ are also shown. The improved convergence of the $H^{-2}$ metric is clear; the benefit of increasing the order of the finite elements themselves is less clear.

Example 7

In this example we study the sensitivity of the norms to small wiggly perturbations. The base shape is the same circle of radius 0.5, and the domain remains $\Omega=[-1,1]^{2}$. The circle is perturbed by scaling its radius by a factor $1+\epsilon\cos(\omega\theta)$; see the illustration in Figure 13. We take 5000 points on the shapes and compute the distances using piecewise linear finite elements and apply Richardson extrapolation to the results from meshsizes $M=80$, $M=160$, and $M=320$ to ensure that finite element discretization errors are negligible. The results for the two norms are shown in the table in Figure 13. As expected from the analysis of nearby straight sections in Section 3.4, for fixed frequency $\omega$ the distances are $\mathcal{O}(\epsilon^{1/2})$ for the $H^{-1}$ metric and $\mathcal{O}(\epsilon)$ for the $H^{-2}$ metric. However, it appears that overall the distances are independent of frequency for $H^{-1}$ and $\mathcal{O}(\omega^{-1/2})$ as $\omega\to\infty$ for $H^{-2}$. Overall, the $H^{-2}$ metric is much less sensitive to small wiggles, especially small high-frequency wiggles, than the $H^{-1}$ metric.

Example 8

In this experiment we consider a family of shapes, the supercircles $x^{r}+y^{r}=(1/2)^{r}$. shapes, together with their representers, are shown in Figure 14 with 512 points on the curves, meshsize $M=80$ and piecewise linear elements, and $r=2^{1},2^{1.5},2^{2},2^{2.5}$. Note that although the different between the curves is apparent, it is much harder to see a difference between the representers. Nevertheless, the norms do see a difference, with the $H^{-1}$ norm decreasing from 4.67 to 3.61 over the four curves shown in the Figure. The $H^{-2}$ norm decreases from 1.86 to 1.26.

Example 9

In this example we study the $H^{-1}$ geometry of the family of supercircles. We take the curves $x^{r}+y^{r}=(1/2)^{r}$ for $0<r<\infty$. This family ranges from a cross (equivalent to a null shape for currents), through an astroid, a circle, to a square. We compute the shapes for 13 exponents $r=2^{j}$, $j=-3,-2.5,\dots,3$. We use 512 points on each curve, a meshsize of $M=10$ and piecewise linear elements. Thus $W$ is a normed vector space of dimension $2\cdot(M+1)^{2}=242$. It is transformed to a standard Euclidean ${\mathbb{R}}^{242}$ using the Cholesky factorization described in Section 3.4. The 12 data points in ${\mathbb{R}}^{242}$ are projected to ${\mathbb{R}}^{2}$ using standard PCA. The results are shown in Figure 15. The bunching up of the points at either end is clear. The induced geometry of the family can be seen in the curvature of the family.

Example 10

In the previous example, the dataset was intrinsically 1-dimensional. We now consider a high-dimensional dataset. We generate 32 random smooth shapes using random Fourier coefficients and then compare them in the $H^{-1}$ metric using currents $\int x^{m}y^{n}\,\textrm{d}x\,\textrm{d}y$ for $0\leq m+n<10$. The Fourier coefficients $\tilde{z}_{k}$ of the shapes have $\tilde{z}_{0}=0$, $\tilde{z}_{1}=0.5$ $\tilde{z}_{-1}=0$ (so that they are all roughly centred and of the same size), with $\tilde{z}_{k}$ for $2\leq k\leq 6$ being independent random normal variables with standard deviation $1/(1+|k|^{3})$. Thus the dataset is 10-dimensional. We then perform an optimization step that computes the planar embedding of the shape currents whose Euclidean distance matrix best approximates the $H^{-1}$ distance matrix, using the best planar subspace of $W$ as computed by PCA as the initial condition for the optimization. The results are shown in Figure 16. The mean distance error of the embedding is 0.06. This gives a pictorial representation of a high-dimensional data set lying within an infinite-dimensional nonlinear shape space. The triangular shapes appear on the edges, sorted by orientation, while the squarish shapes appear near the centre. Obviously similar shapes, such as 5 and 16, and 2 and 4, are placed very close together.

Example 11

In order to compare our approach to the differential signature method, we took the same shapes as in 10 and compared them modulo the special Euclidean group of the plane by computing their Euclidean differential signatures $(\kappa,\kappa_{s})$ (where $\kappa$ is the Euclidean curvature and $s$ is arclength) using the Euclidean-invariant finite differences introduced by Calabi et al. Calabi et al. (1998) and corrected by Boutin (Boutin, 2000, Eq. (6)). The signatures are mapped into $[-1,1]^{2}$ as $(0.8\arctan(\kappa/3),0.8\arctan(\kappa_{s}/150))$; this mapping controls the weighting of the extreme features of the signatures at which $\kappa$ and (especially) $\kappa_{s}$) are large. The currents using moments for order $<10$ are computed as for Figure 16, and the currents embedded in the plane to preserve as best as possible their pairwise distances by performing a least-squares optimisation in this space to minimise the distance between the pairwise distances between the points in the original space and in this 2D version.

There are some interesting differences between the two embeddings, such as the tight clustering of 14, 24, 26, 28 in this Figure, of which only 24 and 26 are close in Figure 16. Looking at the position of shape 10 in both Figures, it is also clear that the differential signature is dominated by the slight kink on the right of the shape (see also shape 13).

Example 12

Our next example (Figure 18) is a family of complicated shapes that have a linear dependence on a parameter, $\phi(t;a)=\phi_{1}(t)+a\phi_{2}(t)$. As the shapes are relatively complicated and the differences between them relatively small, the currents of the shapes (again, moments of order $<10$) embed extremely well in a plane: mean distance error $<10^{-3}$. However, in this embedding in the plane, significant geometry is seen within the family, with speed along the curve decelerating, changing direction, and accelerating. This example and the previous one illustrate that the method of finite element currents copes with complicated immersions.

Example 13

In our final example, we investigate how currents can be used to separate populations of shapes. We generate a set of shapes using six Fourier components as previously, but we choose different values for the 3rd and 4th Fourier components and add relatively small noise to them. The resulting shapes thus come from three different closely-related sets. We generate a small dataset comprising ten examples from each of these sets. These are shown in Figure 19 in the three sets; it can be seen that it is hard to distinguish between the examples.

We then apply PCA to the currents and use PCA to project the shapes into 2D. Figure 20 shows the results of this PCA with a few examples of each shape plotted. The three classes are clearly visible in this embedding.