Dictionary Learning for Two-Dimensional Kendall Shapes
Anna Song Virginie Uhlmann Julien Fageot Michael Unser


Out-of-equilibrium operation of a quantum heat engine:
 The cost of thermal coupling control
    Joachim Ankerhold
    December 30, 2023
============================================================================================




We propose a novel sparse dictionary learning method for planar shapes in the sense of Kendall, namely configurations of landmarks in the plane considered up to similitudes. Our shape dictionary method provides a good trade-off between algorithmic simplicity and faithfulness with respect to the nonlinear geometric structure of Kendall's shape space. Remarkably, it boils down to a classical dictionary learning formulation modified using complex weights. Existing dictionary learning methods extended to nonlinear spaces either map the manifold to a reproducing kernel Hilbert space or to a tangent space. The first approach is unnecessarily heavy in the case of Kendall's shape space and causes the geometrical understanding of shapes to be lost, while the second one induces distortions and theoretical complexity. Our approach does not suffer from these drawbacks. Instead of embedding the shape space into a linear space, we rely on the hyperplane of centered configurations, including pre-shapes from which shapes are defined as rotation orbits. In this linear space, the dictionary atoms are scaled and rotated using complex weights before summation. Furthermore, our formulation is more general than Kendall's original one: it applies to discretely-defined configurations of landmarks as well as continuously-defined interpolating curves. We implemented our algorithm by adapting the method of optimal directions combined to a Cholesky-optimized order recursive matching pursuit. An interesting feature of our shape dictionary is that it produces visually realistic atoms, while guaranteeing reconstruction accuracy. Its efficiency can mostly be attributed to a clear formulation of the framework with complex numbers. We illustrate the strong potential of our approach for the characterization of datasets of shapes up to similitudes and the analysis of patterns in deforming 2D shapes.




	Keywords
	Kendall's shape space; sparse dictionary learning; 2D shape analysis; interpolating curves; splines.





	AMS
	94A12, 65J22, 51M99, 65D18, 65D05





§ INTRODUCTION


Shape analysis is highly relevant to biomedical imaging and computer vision. Among many other applications, it may be deployed to retrieve the main features from a collection of shapes, compare shapes as in morphometrics <cit.>, classify or recognize objects <cit.>, describe the dynamics of a moving object or organism <cit.>, or study the distribution of data in a shape space <cit.> relatively to a mean shape <cit.>.
Depending on the application,
the concept of shape has different meanings <cit.>. Traditionally, shapes are handled as vectors _1,...,_K ∈^d and, for instance, correspond to silhouettes of objects outlined with landmarks. Standard signal analysis tools such as principal component analysis (PCA) can then be used to find the main modes of variation in the dataset of shapes <cit.>.

Alternatively, one can improve results obtained with PCA by using more refined tools, such as sparse dictionary learning  <cit.>. Given a dictionary of representative elements referred to as atoms, a signal is reconstructed from a sparse linear combination of them that minimizes the approximation error. This sparse coding is motivated by the assumption that natural signals are sparse <cit.>. The notion of sparsity has already proven its importance in
an extensive range of problems, from image denoising and signal recovery to recognition and classification <cit.>. Moreover, when the dictionary is learned from the data, significant improvements can be made on the signal reconstruction <cit.>, leading to the so-called sparse dictionary learning approach, classically formulated as

    inf_,∑_k=1^K _k - _k^2 + λSp(_k),

where = (_̣1,...,_̣J) ∈^d × J is the dictionary and = (_1,...,_K) ∈^J × K contains the weights used to reconstruct _k. A bounding constraint on the atoms, _̣j≤ 1, is added to ensure the well-posedness of the minimization problem. Sparsity is promoted by the term λSp(_k) which penalizes nonzero coefficients in the weights, with λ > 0 a parameter. Often, one relies on Sp(_k) = |_k|_0 or Sp(_k) = |_k|_1, which denote the ℓ_0 constraint and the ℓ_1 norm, respectively. When  and  approximately optimize (<ref>), _k is the reconstruction of the original data _k. 

In some applications, the data are pre-processed in order to discard the influence of uninformative features such as the specific position, size, and orientation of the silhouettes. As in Procrustes analysis <cit.>, original silhouettes are scaled, translated, and rotated so as to optimally match a reference silhouette, typically the mean of the dataset.
However, for data with high variability, this approach sometimes fails to produce visually interpretable representative elements (modes or atoms) because, after alignment, silhouettes are still not position-, scale-, and orientation-independent.
To address this, we should handle data as shapes in the sense of Kendall <cit.>. By definition, two geometric objects have the same shape if they are equivalent up to (direct) similitudes (i.e., up to translation, scaling, and rotation).
Shapes are hence considered as nonlinear objects. Our work is motivated by the lack of a simple and efficient sparse dictionary learning method dedicated to Kendall's shape space. Such an analysis would be truly invariant to similitudes, and benefit from both the efficiency of sparse dictionary learning and a valuable shape analysis framework <cit.>.

An additional interesting feature of our approach is that it is suitable for the analysis of shapes defined from discretely-defined configurations of landmarks as well as continuously-defined curves.
In Kendall's original formulation, the shape space is built upon configurations of landmarks in (^d)^N considered up to similitudes <cit.>.
Here, we extend Kendall's shape space to interpolating curves linearly generated by finitely many basis functions.
This extended framework applies in particular to spline curves
generated by piecewise polynomials.
They constitute a common representation for parametric curves <cit.>. Our work is hence related to <cit.>, which relies on the isometry between spline curves and configurations of landmarks to apply dictionary learning after pre-alignment. 




  
Related works

Standard dictionary learning methods cannot be straightforwardly extended to Kendall's shape space because the latter is not a vector space but a Riemannian manifold. The difficulty resides in the nonlinear geometric structure of this space, in which a meaning must be given to linear combinations of atoms. Until recently, most dictionary learning approaches were devoted to data lying in linear spaces. In the past years, however, a few works have focussed on the extension of sparse dictionary learning to nonlinear spaces such as the Grassmann manifold, the manifold of symmetric positive definite (SPD) matrices, Kendall's shape space (as in the two works <cit.> that are closely related to ours), and more general Riemannian manifolds. This can be achieved in two ways: either by mapping the manifold to a Hilbert space, typically a reproducing kernel Hilbert space (RKHS) <cit.>; or by projecting the data onto a tangent space, once at a reference point, or iteratively at multiple points <cit.>. Both approaches flatten the nonlinear space by mapping it to a linear one, so as to make weighted sums of atoms possible.

On the one hand, the mapping of Kendall's space to a RKHS through the kernel trick, as done in <cit.>, follows a classical procedure and enables, after mapping, the use of several established techniques. However, it looses the simple structure of the underlying shape space because the mapping is not explicit. On the other hand, projecting the dataset on one tangent space unfaithfully represents the original distances outside of a neighborhood around the pole <cit.>. In <cit.>, an analysis with multiple projections instead is proposed by the authors. They adapt to Kendall's 3D shape space a sparse dictionary learning method first proposed in <cit.> for general Riemannian manifolds, with the analysis of trajectories of shapes for 3D action recognition as targeted application.


Although the approach of <cit.> has proven its success for action recognition tasks, some of its characteristics are less relevant here. It requires one to fold and unfold the manifold several times on tangent spaces using log__k and exp__k mappings. Most importantly, an affine constraint is added in order to ensure a nontrivial solution, but this requirement modifies the original problem, as pointed out by <cit.>. In addition, the original shapes _k themselves are still needed in their reconstruction, which is undesirable here since the goal is to reconstruct a dataset from the dictionary and weights only. As we do not intend to tackle the same applications, our method should be considered as an alternative for learning shapes. 



  
Overview of our method
In contrast to all previous approaches, the sparse dictionary learning method that we propose is mathematically simple while remaining faithful to the nonlinear structure of the shape space. We avoid kernel methods because they provide no additional geometrical understanding to the problem and are too sophisticated in regard of the simple structure of the shape space. We do not rely on tangent projections either, thereby avoiding both distortions and theoretical complexity.

Instead, the key idea of our approach can be summarized as follows. Kendall's shape space is the quotient of a pre-shape sphere  by the group of rotations. The manifold  is the sphere inside the linear space composed of all centered configurations. We compute all linear combinations in this vector space. Our atoms are pre-shapes in  (i.e., centered and normalized configurations). We use complex numbers in the linear combination so as to scale and rotate them. We then add them together to obtain a configuration whose shape is close to that of the original one.
This proximity is measured through a particular metric in Kendall's shape space. For the classical full Procrustes distance, our method simply leads to the optimization problem

    inf_,  : -.5cm _̣j ∈
     |_k|_0 ≤ N_0 ∑_k=1^K |_k - _k |_^2,

where |·|_ is an ℓ^2 norm specific to the representation (landmarks or interpolating curves), _k are complex weights, and the data _k and atoms _̣j are (complex) pre-shapes. Our work therefore proposes a natural extension of the standard dictionary learning method with a strong theoretical justification.



  
Contributions


	
  * We extend Kendall's original framework to continuously-defined interpolating curves. For this purpose, we put forward the concept of configuration, which indifferently represents a discrete object (configuration of landmarks) or a continuous one (interpolating curve).
	We thus extend and embed the work from <cit.> inside the more general (and appropriate) framework of Kendall's shape analysis.
	
  * Our entire approach of the 2D shape space is formulated in terms of complex numbers and Hermitian inner products. The geometrical interpretation of complex numbers provides a clear understanding of the framework.
	
  * Our main contribution consists in a simple and efficient sparse dictionary learning method, that we call 2D Kendall Shape Dictionary (2DKSD), dedicated to the analysis of 2D shapes in the sense of Kendall. Our approach provides faithful reconstructions with respect to the nonlinear geometric structure of Kendall's shape space, while remaining mathematically light. Using the full Procrustes distance <cit.> to compare the original and reconstructed shapes, 2DKSD boils down to a simple and nearly classical dictionary learning formulation (<ref>) that relies on complex weights instead of real ones. It allows one to scale and rotate the atoms inside the weighted sum individually for each data point, instead of aligning the dataset to a reference shape as a pre-processing step.
	
  * Our implementation of this simple formulation is an adaptation of the algorithm used in the  software <cit.> to the Hermitian framework, and is freely available online.[https://github.com/ansonang3/2DKSD] More precisely, we combine the method of optimal directions <cit.> to an order recursive matching pursuit <cit.> with a Cholesky-based optimization.
	
  * Thanks to the complex setting, our method provides better reconstruction accuracy than approaches relying on the real setting. The atoms of 2DKSD are also visually more realistic and similar to shapes present in the original dataset.




  
Outline of the article

In Section <ref>, we introduce the notion of planar configurations , and describe the action of similitudes over them.
Section <ref> is devoted to Kendall's space of 2D shapes, reformulated for general configurations. We briefly recall the structure of the shape space and its three classical metrics. Then, in Section <ref>, we expose our main contribution: a sparse dictionary learning method dedicated to the analysis of 2D Kendall shapes, which for a well-chosen shape metric leads to a nearly standard formulation with complex weights (<ref>). We expose the implementation of 2DKSD in Section <ref>. Finally, we validate our approach in Section <ref> by experimenting on shapes extracted from real image datasets, before concluding in Section <ref>.














§ CONFIGURATIONS IN THE PLANE AND SIMILITUDES





 §.§ A fundamental example: landmarks (discrete)


Traditionally, the silhouette of a 2D object is sampled by an ordered set of points, called landmarks, that often mark salient features of the boundary. The object itself is then represented by this collection of N points, which we call a configuration of N landmarks, and is an element of ^N denoted by the bold letter . We are now interested in how (direct) similitudes transform configurations. In the sequel, the word similitudes refer to direct similitudes.

It is elementary but important to observe that the action of a similitude in the plane, denoted by (a,b) ∈^* × (where ^* := ∖{0}), has a nice expression using complex numbers

    ∀ z ∈,      (a,b)  ⊙  z = a z + b,

where the multiplication by a applies to z a scaling of modulus |a| and a rotation of angle a, while the addition with b encodes the effect of a translation. We recall that  is the complex argument: for z in , z = |z| e^i(z), with (z) ∈ [0,2 π).
Similitudes form a group (𝒢,∘) whose composition law is given by (a',b')  ∘  (a,b) = (a'a, a'b + b').

Similarly, the action of similitudes on configurations of landmarks also has a natural expression

    ∀∈^N,      (a,b)  ⊙  = a  + b ,̆

where =̆ (1,...,1) is the configuration collapsing to 1. In Figure <ref> (left), the two configurations _1 and _2 are obtained from each other by a similitude. We say that they are equivalent up to similitudes.



 §.§ Interpolating curves (continuous)


We now move to a more involved concept of configuration: interpolating curves, that are used as an alternative to landmarks for representing silhouettes. In the plane, interpolating curves with N degrees of freedom are built by interpolating N complex coefficients [0],...,[N-1] using linearly independent basis functions ϕ_n ∈Ł^2([0,1],), n = 0,...,(N-1). They have the general expression

    ∀ t ∈ [0,1],    r(t) = ∑_n = 0^N-1[n] ϕ_n(t).

If r(0) = r(1), then r is a closed interpolating curve, otherwise it is an open interpolating curve.
The vector ∈^N is called “the” control vector of r (as we show below,  is uniquely defined).


The ϕ_n are often taken to be continuous, as in the case of interpolating spline curves, or simply spline curves, for which the ϕ_n are piecewise-polynomial.
Interpolating curves are considered as continuously-defined objects, in opposition to discretely-defined landmarks. Yet, they are in fact intermediate objects between the continuous and discrete settings. As illustration, we describe the construction of interpolating curves relying on cubic B-spline interpolation in Appendix <ref>.




  
Isometry
Interpolating curves are elements of the space of square-integrable curves
:̋= Ł^2([0,1],).
We endow $̋ with the standard Hermitian inner product

    rs := ∫_0^1r̅  s = ∫_0^1 (r_1 s_1 + r_2 s_2)   + i∫_0^1 (r_1 s_2 - r_2 s_1),      |r|_H = √(rr),

wherer = r_1 + i r_2 ∈$̋, and similarly for s.
Let Γ denote the linear map

    Γ : [                 ^N                  →                                     ↦ ∑_n = 0^N-1[n] ϕ_n ]

and ∈^N × N the Gram matrix of the basis functions ϕ_n, defined by

    [n,m] := ⟨ϕ_n, ϕ_m ⟩_Ł^2([0,1],) = ∫_0^1 ϕ_n(t) ϕ_m(t)  dt.

The matrix  is real-valued, symmetric, positive-definite. 
Let us endow (^N,) with the Hermitian inner product (or shortly, Hermitian product) associated to  : the product of ,∈^N is = ^*, where ^* := ^T refers to the conjugate transpose of . The corresponding norm is then ||_ := √(^*). In Appendix <ref>, we detail how the Hermitian product is related to the real scalar product. As a recall of the conjugate symmetry, please note that ^* = ^*. 

 
	The linear map Γ is an isometry from (^N,) to (,̋  ), in the sense that
	
    ∀ ,∈^N,   Γ()Γ() = .

	Therefore, interpolating curves form an N-dimensional subspace Γ(^N) ⊂$̋ isometric to(^N,).

	We develop the product ∑_n = 0^N-1[n] ϕ_n ∑_m = 0^N-1[m] ϕ_m to obtain 
∑_n,m[n][m] ⟨ϕ_n, ϕ_m ⟩_Ł^2([0,1],), which is exactly ^*.


As a consequence, an interpolating curver ∈Γ(^N)is completely identified to its control vector= Γ^-1(r) ∈^N. 
This result is important for the extension of Kendall's theory to interpolating curves[The question of defining a good notion of curve shape and, in particular, a suitable representation for the curve (such as the square-root velocity (SRV) representation <cit.>), is outside the scope of this article. We refer interested readers to <cit.>.] and ensures that they can be identified to configurations.

  
Action of similitudes over interpolating curves
The action of any similitude(a,b) ∈^* ×on curves in$̋ has once again a nice formulation

    ∀ r ∈,̋     (a,b)  ⊙  r = ar + b 1,

where 1 refers to the constant function ∀ t ∈ [0,1],  t ↦ 1 in $̋.
Suppose that1 ∈Γ(^N), namely, that constant curves are generated by the same basis functionsϕ_n. Under this assumption, it is legitimate to define:̆= Γ^-1(1).
The action of similitudes is induced on^Naccording to

    ∀∈^N,      (a,b)  ⊙  = Γ^-1((a,b)  ⊙ Γ()) = Γ^-1(a r + b 1) = a  + b .̆

This is well defined, becausea r + b 1remains inΓ(^N)for any interpolating curver.



 §.§ A more general concept: configurations


The two previous frameworks - landmarks configurations or interpolating curves - can be encompassed inside one common and more general framework.


	A configuration of dimension N ∈^* in the plane is a complex vector ∈^N.


The notion of configuration plays a major role throughout this article and allows us to handle planar objects, such as curves or landmarks configurations, parameterized by{1,...,N}through the configurations∈^N.

Set an integerN ≥1. Let us consider a Hermitian matrix= ^* ∈^N ×Nthat is positive-definite:∀∈^N ∖{0},^* > 0. Note thatis not necessarily real-valued here. We endow^Nwith the Hermitian inner product := ^*. Working with landmarks in the usual way,refers to the identity matrixId ∈^N ×N. In the case of interpolating curves,is the Gram matrix of the basis functionsϕ_n, assumed to be linearly independent, and to linearly generate1.

Now let$̆ denote a non-zero configuration, that we call the shift configuration, and consider the group action ⊙ : (𝒢,∘) ×^N →^N expressed as (a,b)  ⊙ ↦ a  + b ,̆ where (𝒢,∘) is the group of similitudes. The shift configuration conveys an important geometrical interpretation, and is the vector determining the action of translations on ^N, hence the term “shift”. Together with , they will define the shape space built in next section. In our fundamental examples, we take =̆ (1,...,1) for landmarks, and =̆Γ^-1(1) for interpolating curves. In both cases, $̆ specifies the coefficients of the unit vector or function in the corresponding basis of the representation.

In what follows,$̆ denotes the complex vector line generated by $̆, and()̆^⊥_the complex hyperplane of configurations orthogonal to$̆.


	If ∈$̆, thenis said to be degenerate, and non-degenerate otherwise.

In our examples of landmarks or interpolating curves, degenerate configurations represent objects that collapse to a single point in the plane, and are not interesting in practice. 


	A configuration  is said to be centered if ^̆*  = 0 or, equivalently, ∈ ()̆^⊥_.
	
	Centering consists in orthogonally projecting  onto ()̆^⊥_.

For any configuration, there exists a unique complex numberb ∈, called the center of, so that- b $̆ is orthogonal to $̆. It is given byb = ^̆* /||̆^2_, and- b $̆ is then the centered version of .
In the case of landmarks configurations, the center is the usual arithmetic mean. The operation ↦( - ∑[i]/n) centers . Note that ||̆_ = √(n). Regarding interpolating curves, the center of  coincides with the temporal mean of r, given by
r̅ := ∫_0^1 r(t)  dt ∈ = 1r = ^̆*  = ^̆* /||̆^2_.
Note that, contrarily to landmarks, ||̆_ = |1|_Ł^2([0,1],) = 1.


	A configuration  is called a pre-shape if it is centered and normalized
	
    ^̆*  = 0   and   ||_ = 1.



To visualize the various operations, we encourage the reader to refer to Figure <ref>. We also summarize some notations in Table <ref>.





§ KENDALL'S SPACE OF PLANAR SHAPES



We now describe Kendall's shape space in the planar case <cit.>. When developing this framework, we make a systematic use of the shift configuration $̆ and Hermitian inner productintroduced in the previous section. In particular, interpolating curves can be referred to as configurations. We recall that we consider the action of similitudes written as(a,b)  ⊙  = a + b $̆.



 §.§ Pre-shape sphere and shape space

A shape, in the sense of Kendall, is an equivalence class of non-degenerate configurations considered up to similitudes (see Figure <ref>). It is what remains after discarding redundant geometric information given by the position (center), scaling, and orientation. The shape space is then the quotient of the non-degenerate subspace by the action of similitudes.
We first get rid of translation and scaling, thus obtaining centered and normalized configurations which form the pre-shape sphere (see Figure <ref>). Then, we quotient by the group of rotations
U(1) := {z ∈ |  |z| = 1}
to obtain Kendall's shape space.


	The pre-shape sphere = {∈^N  | ^̆*  = 0, ^*  = 1 } is the set of centered and normalized configurations ∈^N. It is the quotient of ^N by the action of translations ↦ + b $̆,b ∈, and of scalings↦λ  ,λ∈^*.

Since it consists of the unit-norm elements of the complex hyperplane()̆^⊥_orthogonal to$̆,  is a smooth compact real (2N-3)-manifold, as can be seen in

    = {||_ = 1 }  ∩   ()̆^⊥_.


Any non-degenerate configuration is uniquely associated to a pre-shape by centering and then normalizing it (but the converse is not true).

The pre-shape uniquely associated to a non-degenerate configuration ∈^N is the projection of  onto , denoted by

    Π_() = _0/|_0|_,     _0 =  - ^̆*/||̆_^2,̆

where _0 is the centered version of  (see Figure <ref> and Table <ref>).



	The shape space
	
    :=  / U(1)

	is the pre-shape sphere quotiented by the action of the group of rotations ↦e^iθ, θ∈ [0,2π).
	The orbit of a pre-shape ∈ under the action of U(1) is denoted by [] ∈ and is called the shape of . More generally, given a non-degenerate configuration , we define its shape  to be that of the unique pre-shape associated to .


Kendall's shape space  is a compact smooth real (2N-4)-manifold, identified to the complex projective space <cit.>

    ℙ^N-2 = (^N-1∖{0}) / ^* ≃𝕊^2N-3 / U(1).


In this expression, ^N-1∖{0} refers to the set of centered non-degenerate configurations of landmarks. In the literature, a typical approach for quotienting out translations is to discard the last coordinate in ∈^N. For centered configurations of landmarks, one indeed has that [N-1] = -∑_n = 0^N-2[n]. In the general case, ^N-1∖{0} can be replaced by ()̆^⊥_∖{0} and 𝕊^2N-3 by , conserving the identification of  to ℙ^N-2.




 §.§ Distances in the shape space


Three classical distances are usually defined on the shape space: the full (Procrustes), partial (Procrustes), and geodesic distances, denoted as d_F, d_P, and ρ, respectively <cit.>. Their evaluations dist([],[]) enjoy simple and closed-form expressions that involve only the pre-shapes  and . Readers interested in the practical use of these metrics and not in their geometrical definition may directly skip to (<ref>), (<ref>), and (<ref>). For a quick geometrical intuition without mathematical details, we refer readers to Figures <ref> and <ref>.




  
Distances and optimal transformations
In this paragraph,  and  denote pre-shapes by default.



	Let ,∈. We define
	
    d_P([],[]) = min_θ∈ [0,2π) |e^iθ  - |_.

	We say that we optimally rotate  along  when, if unique, the optimal rotation angle θ(,) is applied to  (see Figure <ref>).



	
	Let ,∈. We define
	
    d_F([],[])    = min_a ∈ |a   - |_
       = |P_ - |_,

	where P_ is the orthogonal projection onto the complex vector line generated by  (with respect to ).
	In (<ref>), the optimal alignment factor is given by a(,) = ^*, and its modulus λ(,) is called the optimal scaling factor (see Figure <ref>).
	We say that we optimally align  along  when a(,) is applied to  to obtain P_ = a(,)   = ^*.





 
	Let ,∈ and suppose that ^*≠ 0.
	Then, the optimal rotation angle is given by
	
    θ(,) =  (^*).

	As a consequence, the angle involved in the optimal alignment factor is the same as the optimal angle itself: (a(,)) = θ(,).
	Otherwise, when ^* = 0, the distances are maximal and equal to d_F([],[]) = 1 and d_P([],[]) = √(2).


	If ^* = 0, then we use the orthogonality of  and  to conclude. Otherwise, since
	
    _θ∈ [0,2π)  |e^iθ· - |^2_ = _θ∈ [0,2π)≤ |e^iθ·|e^iθ·,

	θ is optimal when e^iθ· = e^-iθ∈^+; in other words, when θ =.



Note that the orthogonality condition ^* = 0 means that  and  are decorrelated: whatever rotation we apply on , the distance |e^iθ· - |_ remains maximal. In the real setting, this corresponds to the case where  is orthogonal to all rotations of  with respect to the real inner product (see Appendix <ref>).



	
Let ,∈. The distances have closed-form expressions given by

    d_P([],[])    = |^*/|^*| - |_ = √(2 - 2 |^*|)∈ [0,√(2)], 
    
    d_F([],[])    = |(^*)   -  |_ = √(1 - |^*|^2)∈ [0,1].



In Figure <ref>, we illustrate the fact that the distances d_P and d_F measure the shortest reachable norm | - |_, if  is the image of  after any rotation, or scaling and rotation.



	
	Let ,∈. The Riemannian distance in the shape space (see Figure <ref>), also called geodesic distance and denoted by ρ([],[]), is equal to
	
    ρ([],[]) = arccos |^* | = arccosλ(,)   ∈ [0,π/2].








	The three distances d_F, d_P, and ρ are related to each other as
	
    d_P^2    = 2 - 2 √(1 - d_F^2), 
    
    	d_F    = sinρ, 
    
    	d_P    = 2 sin (ρ/2).




We found it worth presenting the Riemannian distance ρ alongside d_F and d_P, although our 2DKSD method does not require the use of the Riemannian structure. In fact, the geometric structure of  is entirely inherited from that of : this justifies using the pre-shape sphere, whose geometry is well understood, to think about the shape space itself. In particular, the Riemannian structure of  is a quotient structure of , as discussed in Appendix <ref>. As a consequence, tangent spaces, geodesics, and exponentials benefit from explicit expressions on . For readers interested in the Riemannian structure of the shape space (not restricted to the planar case), we refer to <cit.>, where the usual definition relying on landmarks is used. Examples of more advanced references about Riemannian geometry, in particular that of complex projective spaces, are <cit.>.




§ SHAPE DICTIONARY


The simplest dictionary learning approach on a dataset of configurations consists in considering them as elements of a point cloud in (^N)^K ≃ (^2N)^K, and then applying an existing dictionary learning method. However, as discussed in the introduction, a pre-processing of the data can be necessary in order to attenuate differences of position, scaling, or orientation across the dataset. These geometric dissimilarities indeed lead to undesirable artefacts after learning, such as smoothed and unfaithful reconstructions or distorted and redundant atoms (in the sense that they share roughly the same shape but are associated to e.g. different orientations).
An example of pre-processing is proposed in <cit.>, where the data are first optimally aligned
along a mean configuration _mean, with _mean∈ being the Fréchet mean of _1,...,_K∈^K with respect to the distance d_F (see Appendix <ref>).

A limitation of this kind of “align-first” approach is that the alignment step depends on the whole dataset, and needs to be started anew for each new data input. In Section <ref>, we compare the results obtained by aligning first to those produced by our algorithm. Instead, we prefer to directly work in the shape space with suitable shape metrics, without having to pre-process the data. 

Our 2DKSD method takes advantage of complex weights to rotate and scale the atoms before summing them to reconstruct each original shape. Hence, instead of pre-aligning the data, we rotate and scale inside the weighted sums _k, individually for each data shape. Also, we show that, for a good choice of the error metric that serves to compare the original and reconstructed shapes, the 2DKSD boils down to a nearly classical and very simple dictionary learning formulation. Weights are however complex vectors, and the dataset and dictionary atoms must satisfy some (light) constraints imposed by Kendall's framework. This a priori unexpected result, shown in Proposition <ref>, is interesting both mathematically and numerically.



  
Statement of the Shape Dictionary
We suppose that elements of the dataset are pre-shapes _1,...,_K ∈ (centered and normalized). Let dist denote any distance on , such as one of the three distances d_F,d_P,ρ.

The 2DKSD general formulation of the problem to be solved writes

    inf_, : _̣j ∈,
    _k ∈
     |_k|_0 ≤ N_0    ∑_k=1^K dist([_k], [_k])^2.


The dictionary corresponds to = (_̣1,...,_̣J) ∈^N × J, whose atoms _̣j ∈ are pre-shapes. Importantly, the weights _1,...,_K ∈^J are complex numbers
which apply scalings and rotations to the atoms _̣j before summing them to _k.
This linear combination is sparse, since the (hard) sparsity constraint |_k|_0 ≤ N_0 enforces that at most N_0 coefficients are non-zero. We also impose _k and _k to be pre-shapes before comparing their shapes. If  and  are (approximately) optimizers of (<ref>), the original shape [_k] is reconstructed[When discussing the sparse coding step, we shall show that, for dist = d_F, one can compute _k so that the original pre-shape itself is reconstructed as _k.] as [_k].

The unit-norm constraints imposed on _̣j and _k are in fact nonessential. Up to a rescaling of the weights _k, we can more generally consider _̣j as a centered configuration of norm smaller than 1, thus allowing possibly collapsing atoms _̣j = 0. This does not change the problem since this corresponds to a case where _̣j would be ignored in the linear combinations of the original optimization problem (<ref>). Hence, we could replace _̣j ∈ by _̣j ⊥_;̆ |_̣j|_≤ 1. We have decided instead to keep atoms as pre-shapes to have a uniform scaling across the dictionary.

Also, _k only needs to be a non-degenerate configuration in ^N ∖$̆ for the shape_kto be defined. Since each_̣jis already centered,_kis necessarily centered, and so we only need to suppose that_k ≠0does not collapse. We also note thatΠ_(_k) = _k/|_k|_is then simply the normalized version of_k. Hence, the infimum in (<ref>) is equal to

    inf_, : _̣j ∈
    _k ∈
     |_k|_0 ≤ N_0    ∑_k=1^K dist([_k], [_k])^2 = inf_, : _̣j ∈
    _k ≠0
     |_k|_0 ≤ N_0    ∑_k=1^K dist([_k], [_k/|_k|_] )^2.

Using the full Procrustes distance dist = d_F as minimization criterion, (<ref>) leads to the 2DKSD simple formulation
    inf_ : _̣j ∈∑_k=1^K min__k : |_k|_0 ≤ N_0 |_k - _k |_^2,

where we recall that _k are pre-shapes. This formulation is equivalent to the one in (<ref>).

We express the minimization problem (<ref>) as

    inf_ : {_̣j ⊥_ |_̣j|_≤ 1. ∑_k=1^K inf__k : {_k ≠0
     |_k|_0 ≤ N_0. d_F([_k], [_k/|_k|_] )^2

since, after fixing , the minimization breaks into K independent elementary terms that correspond to a sparse coding of the shapes [_k]. Each term has the form

    inf_ : {≠0
     ||_0 ≤ N_0. d_F ([], [/||_]  )^2 = inf_ : {≠0
     ||_0 ≤ N_0. | - P_|^2_ = inf_ : ||_0 ≤ N_0 | - P_|^2_,

where we successively applied (<ref>), used that the vector line generated by /||_ and  is the same, and dropped[The case = 0 does not induce a lower infimum value because | - P_|^2_ = ||^2_ is then maximal.] the inequality ≠0.
Then, we apply the results of Lemma <ref> and find that

    min_ : {≠0
     ||_0 ≤ N_0. d_F ([], [/||_]  )^2     = min_ : {≠0
     ||_0 ≤ N_0. | - P_|^2_ = min_ : ||_0 ≤ N_0 | - P_|^2_
        = min_||_0 ≤ N_0 | - |^2_,
 where the three first terms share a common minimizer. Therefore, a new formulation of the original problem is given by (<ref>) and, equivalently, by (<ref>).
  

§ DICTIONARY UPDATE AND SPARSE CODING


We numerically solved the 2DKSD problem (<ref>) involving the full distanced_F. Our code is fully available online.

As summarized in Algorithms <ref> and <ref>, the implementation is directly adapted from the algorithm used in the  software <cit.> to a complex Hermitian framework. Recall that= (_̣1,...,_̣J)is the dictionary and= (_1,...,_K)contains the weighting complex vectors.
To minimize the corresponding loss functional

    E(, ) = ∑_k = 1^K |_k - _k|_^2,

we use a classical procedure. After initializing the dictionary with random elements of the dataset, the algorithm alternates between a sparse coding step in whichminimizes (<ref>),
and whereis fixed <cit.>, and a dictionary update step whereis fixed.



 §.§ Dictionary update with MOD

To update the dictionary, it is natural to think about performing gradient descent overfor updating the dictionary. Let us denote by= (_1,...,_K) ∈^N ×Kthe matrix whose columns contain the dataset. The gradient can be conveniently explicited as_E(,) = 2 (^* - ^* ).Here, however, we prefer to rely on the Method of Optimal Directions (MOD) introduced by <cit.>, as done in the  software <cit.>. Given a fixed value of, the dictionary is updated so that it solvesmin_E(,), and its columns are then projected back onto. It is known that, in the landmarks case= Id, a particular solution to the least-squares problemmin_∈^N ×J E(,)is

    = ^+ = Σ^+^*,

where^+ = Σ^+ ^*is the pseudo-inverse of, and the SVD decomposition ofis= Σ^*. The pseudo-inverseΣ^+corresponds to the matrixΣwhere non-zero diagonal elements are replaced by their multiplicative inverse.
We assert in Lemma <ref> of Appendix <ref> that, interestingly, when≠Id, the update is the same and does not depend on.

Let us remark that the solution (<ref>) cancels the gradient (one can see it by checking that^+ ^* = ^*), and also that the update does not involve any previous value assigned to. When(^*)is invertible[This happens whenever the rank of  is equal to its number of lines or, equivalently, when the columns of ^* are independent.], we have^+ = ^* (^*)^-1and the solution is then unique and equal to= ^* (^*)^-1.

After this first operation, the non-zero columns of the dictionary are projected back onto the pre-shape sphere by applying the normalization∀ j = 1,...,J,   _̣j ←_̣j / |_̣j|_.This is sufficient as (<ref>) already centers the atoms, since the data themselves are centered. We replace columns which have become zero, as well as under-utilized ones, by pre-shapes randomly picked in the original dataset.



 §.§ Sparse coding with ORMP

The sparse coding problem for a data pointis formulated as

    min_∈^J 
     ||_0 ≤ N_0 | - |^2_,

whereis known. Equivalently, the non-zero coefficients of the solutioncorrespond to a set of indicesÎthat solves

    min_I ⊂{1,...,J}
     |I| ≤ N_0 | - P_{_̣j}_j ∈ I()|^2,

where{_̣j}_j ∈Idenotes the complex vector space generated by the columns indexed byI.

Let us suppose that a solutionhas been found. Then, the shape[]is approximately reconstructed as[/||_], whenever ≠0. Note that, in fact, the original pre-shapeitself (with preserved orientation) can also be approximated as/| |_: there is no need to rotatein order to bring it close to the original pre-shape. As stated in the proof of Lemma <ref>, this is because the solutionsatisfiesP_= , implying thatcorresponds to the optimal alignment of the pre-shape/||_along(see Definition <ref>). Also,/||_is optimally rotated along(see Proposition <ref>).

  
Order Recursive Matching Pursuit
A possible way to approximate the solution of (<ref>) is to rely on a complex version of Orthogonal Matching Pursuit (OMP) <cit.> adapted to the Hermitian framework(^n,). Instead, we prefer to use Order Recursive Matching Pursuit (ORMP), with a Cholesky-based optimisation <cit.>, as it was empirically found to be more efficient than an adapted complex OMP. The ORMP algorithm is explained in <cit.> and a (non-optimized) implementation is presented in Algorithm <ref>. ORMP starts with an empty set of indicesI[0] = ∅, a remainder^(0) = , and a family of vectors{_1^(0),...,_J^(0)} = {_̣1,...,_̣J}. At each step, the algorithm increasesI[l-1]by a new elementj[l]maximizing|_j^(l-1) * ^(l-1) ||_j^(l-1)|_(wherever defined), in such a manner thatI[l] = I[l-1] ∪{j[l]}minimizes (<ref>) under the nesting constraintI[l-1] ⊊Iand|I[l]| = l.
Then, it updates the vectors{_1^(l),...,_J^(l)}in a way that is similar to the Gram-Schmidt process[The difference being that the order is imposed by j[1],...,j[N_0], so that the final orthonormal family is {_j[1]^(0),...,_j[N_0]^(N_0-1)}.].
It stops either whenN_0indices have been found, or whenis spanned byl_maxindependent atoms and a new indexj[l_max + 1]does not contribute in minimizing (<ref>) further. Finally, the output weightis solution toP_{_̣j}_j ∈I[l_max]() = _I[l_max] _I[l_max]. Its non-zero coefficients are given by

    _I[l_max] = (_I[l_max]^* _I[l_max])^-1_I[l_max]^*  .


  
Optimizing ORMP
Thanks to the so-called Cholesky optimisation of ORMP <cit.>, it is in fact possible to speed up computations considerably. The Gram matrix^* is pre-computed before coding the data(and re-used for other data). At each step1 ≤l ≤N_0, only the elements_j[l]^(l-1) *_̣j /|_j[l]^(l-1)|_
= _j[l]^(l-1) *_j^(l-1)/|_j[l]^(l-1)|_,
     
|_j^(l) |_^2,
    _̣j^* ^(l-1),
  and  _j^(l-1) *^(l-1)/|_j^(l-1)|_
= _̣j^* ^(l-1)/|_j^(l-1)|_are updated, forj = 1,...,J.
The coefficients of the decomposition of^(l) = _j[l]^(l-1)/|_j[l]^(l-1)|_in the linearly independent family{_̣j[1],..., _̣j[l]}are also iteratively computed. If these coefficients are stacked into a column vector^(l) ∈^N_0where the lastN_0 - lcoefficients are zero, the matrix𝐔 = (^(1),...,^(N_0)) ∈^N_0×N_0is upper triangular with strictly positive real scalars on the diagonal. Furthermore, it is the inverse conjugate of the Cholesky factor𝐋involved in the factorisation_I[N_0]^* _I[N_0] = (𝐔^*)^-1 𝐔^-1 = 𝐋 𝐋^*. The final weightis then given by_I[N_0] = 𝐔 [_j[1]^(0) *^(0)/|_j[1]^(0)|_,...,_j[N_0]^(N_0-1) *^(N_0-1)/|_j[N_0]^(N_0-1)|_]^T.All computations can be inferred from the updates in Algorithm <ref> and are not detailed here. They can be transparently investigated in the available source code.



§ COMPUTATIONAL RESULTS


We ran our algorithm on5datasets that vary in size, type and shape.[Experimental results and datasets from this section can be found at https://github.com/ansonang3/2DKSD.]
  – Dataset 1[Images were extracted from videos available from the C. elegans behavioral database <cit.>.] consists of K > 5500 skeletons of Caenorhabditis elegans nematodes delineated as discretely-defined configurations of N = 20 landmarks. The dataset results from the concatenation of videos featuring freely crawling nematodes aged between 2 and 18 days and presenting various locomotion behaviors. The dataset was reconstructed using N_0 = 5 atoms out of the J = 10 atoms learned by 2DKSD.

  – Dataset 2[Data courtesy of <cit.>.] contains K > 6300 skeletons of 4Caenorhabditis elegans nematodes delineated as continuously-defined open Hermite-spline curves <cit.> with N = 12 degrees of freedom. The animals lied in a shared container and were constrained by the lack of space, leading to looping or wavy shapes. The dataset was reconstructed using N_0 = 5 atoms out of the J = 10 atoms learned by 2DKSD.

  – Dataset 3[Data courtesy of <cit.>.] features K = 40 hands of 4 different people, as outlined by N = 56 landmarks. The dataset was reconstructed using N_0 = 3 atoms out of the J = 10 atoms learned by 2DKSD.

  – Dataset 4 and 4b[Data extracted from the binary images of the Kaggle leaf dataset https://www.kaggle.com/c/leaf-classification.] are constituted of K > 1500 leaves outlines, and come in the form of configurations of N = 200 landmarks and cubic B-splines <cit.> with N = 40 degrees of freedom, respectively. For most shapes, the first landmark either marks the stem or the leaf tip. Both datasets were (independently) reconstructed using N_0 = 5 atoms out of the J = 20 atoms learned by 2DKSD.

  – Dataset 5 is a synthetic dataset generated by deforming original silhouettes extracted from the binary images of the MPEG-7 database[Data courtesy of <cit.>. See http://www.dabi.temple.edu/ shape/MPEG7/dataset.html]. For each original pre-shaped configuration _or = (_or^1,...,_or^N) ∈^N, N random Gaussian planar vectors were generated and smoothed with the kernel K = (|_or^i - _or^j|^2)_i,j of the squared distances between landmarks. These deformations were then added to _or. The process was repeated 10 times, resulting in a dataset of K = 14000 shapes featuring 70 classes of objects, each containing 200 configurations of N = 200 landmarks. The dataset was reconstructed using N_0 = 30 atoms out of the J = 230 atoms learned by 2DKSD.

  
Comparison to the “align-first” method 

The problem, that we call “align-first”, relating most to the one we are solving, is stated as

    inf_', ' : -.5cm '̣_j ∈
     |'_k|_0 ≤ N_0 ∑_k=1^K |_k - ' '_k |_^2,     where '_k ∈^J,

where' = ('_1,...,'_K), and_1,...,_K ∈are pre-shapes optimally rotated along a reference pre-shape. Typically, we chose to rotate them along a mean pre-shape_mean(see Appendix <ref>).
This standard problem is essentially the same as in (<ref>), except that the original dataset is pre-processed and, more notably, that the weights'_kare real vectors.

In practice, we took advantage of the MOD combined with the real version of the Cholesky-optimised ORMP provided by the  toolbox <cit.>. This is an especially appropriate comparison, since our algorithm is a direct adaptation of the latter to the Hermitian framework. Thanks to an isometric transformation through√(), (<ref>) could be reformulated as a standardℓ^2problem.

  
Comparison to complex PCA
As a valuable and complementary comparison, we also considered the complex version of principal component analysis. It also relies on complex weights, but atoms are chosen in a specific way, that we recall here in our setting. Suppose that the average configuration has been discarded in each element of the dataset, and let us still use_1,...,_Kto denote this new dataset (the average configuration coordinates are simply the coefficient-by-coefficient arithmetic means). We find the first complex PCA mode by solving

    _|_1|_≤ 1∑_k = 1^K |_k - (_1^* _k) _1|_^2 = _|_1|_≤ 1∑_k = 1^K |_k|_^2 - |_1^* _k|^2 = _|_1|_≤ 1∑_k = 1^K |_1^* _k|^2.

Let us set= (_1,...,_K)as before. The matrix^* is self-adjoint with respect to the Hermitian product, because∀,∈^N, ^*    = ^*   . By the Spectral Theorem, this implies that there exists a-orthonormal basis∈^N ×Nwith^* = Idin which^* is expressed as a diagonal matrix:^*  = ^*  = ^-1.Since^* is positive semi-definite with respect to, we know that= diag(λ_1,...,λ_N)has (real) non-negative values that can be ordered asλ_1 ≥... ≥λ_N ≥0.
A solution of (<ref>) then corresponds to a unit-norm eigenvector of^* corresponding toλ_1.
The second complex PCA mode_2, is the solution to a problem analogous to (<ref>), except that we add the constraint_2 ⊥__1. It is given by a unit-norm eigenvector corresponding toλ_2. Further modes are iteratively found with the constraint to be-orthogonal to previous modes.

The dictionary is then the collection of the modes(_1,...,_J)and the bestN_0-term approximation of_kis∑_j = 1^N_0 (_j^* _k) _j. The corresponding loss is then equal to∑_k = 1^K |_k|_^2 - λ_1 - ... - λ_N_0. As in the previous paragraph, we employ√()to transform the problem into a standard form suitable for computations.



 §.§ Results

In Figures <ref>, <ref>, <ref>, <ref>, <ref> and <ref>, we show the dictionary learned by 2DKSD and the align-first method with specific parameters for each dataset (Table <ref>), alongside with four examples of reconstructed data superimposed over the original one. In each case, the reconstruction errors|_k - _k|_ = |_k - _k|_/|_k|_are expressed in percentage, and can be understood as an absolute error multiplied by100, or as a proportion of|_k| = 1as well. Note that, as mentioned in Section <ref>, the 2DKSD ensures that these errors are equal tod_F([_k],[_k]), which is not the case for the align-first method. We also indicate the root mean square reconstruction errors (RMSE), defined as√( E/K) = √(1/K∑_k=1^K |_k - _k |_^2 )  and  √( E'/K) = √(1/K∑_k=1^K |_k - ' '_k |_^2 ),respectively, and also expressed in percentage. We summarize them in Table <ref>, where we specify the parameters used in the examples, as well as a typical run-time[All experiments were run on a standard laptop with an Intel® Core™ i5-7200U CPU running 4 cores at 2.50GHz, with 7,7 Gb of RAM.]
forT = 30iterations in Algorithm <ref>, which is sufficient for the loss to converge in all datasets, and including the duration of the final sparse coding of the whole dataset. To reduce ORMP run-time for larger datasets (such as Datasets 1 and 2, marked with a star * in the following table), we applied Algorithm <ref> in parallel and independently on each data point_k. To learn overcomplete dictionaries on very big datasets (such as Dataset 5, marked with two stars ** in the table), we chose to update the weights_kstochastically and in parallel by batches of4096, after initializing them to zero. Obviously, additional efforts could be done to optimize the implementation and reduce run-time on big datasets, but this falls out of the scope of the present work, and is kept for the future.

In Figure <ref>, we illustrate in logarithmic scale the RMSE obtained with complex PCA v.s. (real) align-first v.s. 2DKSD, after learning atoms from Datasets 1 to 4b. For varying values ofN_0, we compare the RMSE of a reconstruction relying on the firstN_0PCA modes, and the RMSEs of an align-first or 2DKSD dictionary learned with the parameters(N_0,J).J ≥N_0is fixed to the same values used in the examples shown in the figures, and is such thatJ < N. This comparison is relevant for two reasons. First, the dictionary of PCA atoms is a natural candidate when learning an undercomplete dictionary containing less atoms than the dimension of the configuration space. Second, while complex PCA relies on complex weights but determines the atoms in a specific way, the align-first method is less constrained in the choice of atoms but relies on real weights. Our 2DKSD method is therefore a good compromise between them. We have also explored theJ > Ncase, which leads to an overcomplete dictionary (Dataset 5, Figure <ref>).

  
Comparison with align-first and complex PCA
We argue that the align-first method provides less satisfying results than our shape dictionary, both mathematically and numerically. Mathematically, the optimized loss functional (<ref>), denoted byE', is necessarily larger than our lossE(<ref>) since weights are constrained to lie on the real line. Intuitively, they can scale but not rotate the atoms before summing them. In contrast, the complex weights in our method scale and rotate the atoms, resulting in a smaller lossE. Numerically, this leads to a visible difference in the RMSEs, as seen in Figure <ref>. In particular, the loss obtained with our algorithmEcan be a third or a fourth of the one resulting from align-firstE', as in Figures <ref> and <ref>.
Besides, the dictionary computed by align-first often contains distorted, irregular or unrealistic atoms, which reflect the effort required to fit the data, in spite of the pre-alignment. When analyzing datasets of shapes, it is sometimes preferable to obtain realistic atoms and visually accurate reconstructions instead of perfect reconstructions but unrealistic atoms. Visually realistic atoms indeed offer a way to hypothesize on the nature of the variability of the dataset.

Figure <ref> demonstrates the benefits of using complex numbers in the linear combinations. For increasing values ofN_0, the align-first RMSE indeed becomes significantly larger even when compared to that of complex PCA. As another illustration of their relevance, Figure <ref> features two different weighted sums of the same three atoms that allow reconstructing very dissimilar shapes. It is not surprising for the first hand to be correctly reconstructed as, at first order, its shape is similar to that of the first atom, which is weighted by the coefficient of largest magnitude. The remainder is compensated by the other two atoms. For the second hand, however, it is more surprising to see that hands with open fingers, scaled, rotated, and then summed together, manage to produce a hand with closed fingers.

Last, Figure <ref> shows that 2DKSD is more appropriate than complex PCA forN_0 < J. When learning a KSD dictionary that uses all theN_0 = J < Natoms to reconstruct data, the RMSE is close to that of complex PCA. In fact, the minimal RMSEs should mathematically be equal whenever the average configuration of the dataset corresponds to the zero configuration, or if it is not discarded from the dataset (see the details of complex PCA discussed above). Indeed, both solutions find theJ-dimensional subspace that is the closest to the data points, i.e., the subspace spanned by theJPCA modes.

  
Comments on the examples comparing 2DKSD and align-first
As is seen in Figures <ref> and <ref>, nematode shapes seem to be efficiently reconstructed using not more than5atoms. For both datasets 1 and 2, the reconstructions are visually satisfying and numerically accurate. The results for Dataset 2 are remarkable in that, although the original shapes do not have the characteristic smoothness of freely moving worms, the algorithm has less difficulty in the reconstruction than align-first. Moreover, it produces realistic atoms, whose curve parameterization is well-balanced, in contrast to those of align-first where consecutive control points are not regularly distributed along the curve.

In Figure <ref> reporting the results for Dataset 3, both methods result in very similar losses, but the atoms are less realistic with align-first. They exhibit self-intersecting and stretched fingers. In this case, it is an informative feature as it indicates to which extent fingers “tend” to be closed. Our algorithm is then an alternative to the align-first method.

The results in Figures <ref> and <ref> on the leaf datasets 4 and 4b are similar regardless of the representation (landmarks or B-splines) used. Our lossEis less than84%of the comparison lossE', and we obtain a dictionary with realistic leaves, while align-first contains twisted shapes.

Finally, we observe that the overcomplete dictionaries learned by both methods on Dataset 5, using20of their atoms shown in Figure <ref>, look fairly similar. Thus, 2DKSD offers no particular advantage over align-first regarding the visual realisticness of the atoms in that case, but leads to a significantly smaller RMSE. It also reproduces better high-frequency features of the silhouettes, as in the cow and lizard shapes.




§ CONCLUSION


 Our 2D Kendall Shape Dictionary approach is a natural adaptation of usual learning techniques to Kendall's nonlinear manifold and does not require sophisticated operations. For the classical full Procrustes metric, it simplifies to a nearly standard formulation in which complex weights are used when combining the atoms of the dictionary. As our main contribution, this formulation has the double property to enjoy both a strong theoretical justification and a simple algorithmic framework.
 We have demonstrated the positive impact of using a complex setting: datasets do not require to be aligned along a reference mean because atoms are freely scaled and rotated before being summed to reconstruct shapes independently. This flexibility increases reconstruction accuracy and allows dictionary atoms to remain visually realistic. We have also extended Kendall's space of planar shapes, initially defined for discrete configurations of landmarks, to continuously-defined interpolating curves, by introducing the general notion of configuration. Hopefully, our method is a promising tool for characterizing complex phenotypes from biological images.
 



§ ACKNOWLEDGMENTS

A.S. is especially thankful to Jean Feydy for drawing her attention to reference <cit.>. She is grateful for the support from École Normale Supérieure and École polytechnique fédérale de Lausanne. V.U., J.F., and M.U. were supported by the Swiss National Science Foundation under Grant . V.U. was also partly supported by EMBL core fundings. J.F. was also partly supported by the Swiss National Science Foundation with grant agreement .



§ APPENDIX


 §.§ Hermitian products hold twice more information than scalar products

In the sequel, we make the identification^N ≃^2Nthrough the mapping↦z = ((),())^T, where() = (([0]),...,([N-1]))and likewise for the imaginary part. We use bold and italic to distinguish the complex and real counterparts of the vectors. If(^N,)is endowed with a Hermitian inner product(,) ↦( |  )_:= ^*associated to the matrix∈^N ×N, then one can define on^2Nthe scalar product⟨·, ·⟩_canonically associated to it as

    ⟨ z, w ⟩_ := ( ( | )_),

defining the same norm⟨z, z ⟩_= ( |  )_.Due to the sesquilinearity of the Hermitian product, we have that (i  |  )_= -i ( |  )_, so that ( (i  |  )_) = ( ( |  )_) , and then

    ( | )_ =  ⟨ z, w ⟩_ + i⟨ R_π/2 ⊙  z, w ⟩_,

where byR_θ⊙z = cosθ  ((),())^T + sinθ  (-(),())^Twe denote the image (in the real setting) ofby a rotation byθ, namely, the real counterpart ofe^i θ . Here,R_π/2  ⊙  z = (-(),())^T.

As a consequence,( |  )_= 0if and only if⟨z, w ⟩_= 0and⟨R_π/2  ⊙  z, w ⟩_= 0. In fact, a stronger property can be deduced:∀θ∈[0,2π),⟨R_θ  ⊙  z, w ⟩_= ⟨z, R_θ  ⊙  w ⟩_= 0.
This is becauseR_θ  ⊙  z = cosθ z + sinθ R_π/2  ⊙  z, and we conclude by bilinearity of the scalar product, and then by exchanging the roles ofzandw.
In other words, Hermitian orthogonality can be understood as a real orthogonality between any rotated image of  and  (and conversely).


Suppose that = 𝐀 + i𝐁 with 𝐀,𝐁∈^N × N.
The scalar product is computed as ⟨ z, w ⟩_ = z^T Ψ w,
with Ψ = [  𝐀 -𝐁;  𝐁  𝐀 ].


 §.§ Two lemmas

In this section, we establish two lemmas. The first one is used in the proof of Proposition <ref>, and the second one states that the dictionary update (<ref>) does not involve the matrix.



	Let ∈^J be a minimizer of min_||_0 ≤ N_0 | - |_.
	If ≠0, then any term in (<ref>) is also minimized by .
	Otherwise,  is orthogonal to the columns of  and the distances in (<ref>) are all maximized to 1. In that case, any  satisfying ≠0 is a minimizer of any term in (<ref>).

	The quantity | - |_ reaches a minimum in { |  ||_0 ≤ N_0 }. It corresponds to the minimal distance to  of a subspace generated by at most N_0 atoms _̣j (see previous discussions about sparse coding). Let  denote the corresponding minimizer. It satisfies min_||_0 ≤ N_0 | - |_ = | - |_.
	By definition of a projection, | - | ≥ | - P_|. Moreover, if λ∈ is such that P_ = λ, then | - P_| = | − (λ)|, where |λ|_0 ≤ N_0. The last term is then greater or equal to min_||_0 ≤ N_0 | - |_, hence it comes that
	
    | - |_ = | - P_|_.

	Using the same arguments, one easily shows that this also corresponds to min_||_0 ≤ N_0 | - P_|_.
	As a useful remark, notice that the equality (<ref>) implies that P_ =.
	
	Coming back to the original equalities (<ref>), if ≠0, then all terms are minimized by . Otherwise, if = 0, then  is orthogonal to the vectors _̣j, and  is not correlated to any shape generated by the _̣j, inducing a maximal distance d_F. In particular, any  such that ||_0 ≤ N_0 and ≠0, can be taken as a minimizer of any term in (<ref>).
E being defined in (<ref>), a solution to min_∈^n × J E(,), for a fixed value of , is given by (<ref>), and its expression is independent from the matrix .

We re-write the problem into a standard ℓ^2 form:

    min_ E(,) = min_∑_k = 1^K |α_k - _k|_^2     = min_∑_k |√()α_k - √()_k|^2 
        = min_√()( - ) _Fro^2 
        = min_ (^* ^* - ^*)√()_Fro^2 
        = min_𝐇∑_i=1^n |^* 𝐡_i - 𝐠_i|^2,
𝐡_i and 𝐠_i denoting the columns of 𝐇 = ^* √() and 𝐆 = ^* √(). It is known that the elementary problems are solved by ĥ_i = (^*)^+ 𝐠_i. This gives Ĥ^* = 𝐆^* ^+ and finally = ^+, after inverting with (√())^-1.


 §.§ Mean shape


The mean shape of a dataset can be defined as the Fréchet mean of the points[_k]scattered on Kendall's manifold <cit.>, with respect to one of the three distancesd_F,d_P, orρ. It is the unique global minimizer_meanof∑_k=1^K dist([_mean],[_k])^2,when it exists. The Fréchet mean with respect to the full distancedist = d_F(see Figure <ref>) can be found as the shape of the eigenvector associated to the greatest eigenvalue of the operator <cit.>∑_k P__k,where_k ∈are pre-shapes, and whereP__k = _k _k^* is the orthogonal projector onto the complex vector line generated by_krelatively to the Hermitian product. To prove it, let us consider the problem in the pre-shape sphere. We obtain:

    _∈∑ | - P__k|_^2 = _∈∑ |P__k|^2_ = _∈^* ∑_k P__k,

where we used that|- P__k |_^2 = 1 - |P__k |^2_.
A similar proof for landmarks (= Id) can be found in <cit.>.

  
Link with the mean-shape curve of <cit.>
In this article related to ours, a mean-shape curve of a family of curves{r_1,...,r_K}inH := Ł^2([0,1],^2)endowed with the usual norm is defined as an optimal curve

    r_mean∈_|r|_H = 1 
    r̅ = 0 ∑_k=1^K |r - P_k r|^2_H = _|r|_H = 1 
    r̅ = 0 ∑_k=1^K |P_k r|^2_H,

whereP_kis the similarity projector, i.e., the orthogonal projection onto the subspaceS_r_kof dimension4associated tor_kcontaining all the images up to similitude transforms ofr_k. It writes


    S_r    = {λ[   cosθ - sinθ;   sinθ   cosθ ] r + [ α; β ] | λ∈, θ∈ [0,2π), α, β∈}
        =
    {[ r^x; r^y ], 
    [ -r^y;  r^x ],
    [ 1; 0 ],
    [ 0; 1 ]}.

It can be shown thatr_meanbelongs to the eigenspace associated to the second greatest eigenvalue of∑_k P_kas

    ∑_k=1^K P_k r_mean = λ_2 r_mean,      |r_mean| = 1

and, as a consequence, the mean-shape curve of interpolating curves is also an interpolating curve (as soon as constants are themselves interpolating curves), sincer_mean ∈Im(∑P_k). Aligning the dataset to the mean-shape curve then consists in takingr_k := P_k r_mean.In fact, in the complex setting, ifr_1,...,r_K ∈Ł^2([0,1],)are centered and normalized, the orbit of the mean-shape interpolating curve[r_mean]up to rotations is then exactly the Fréchet mean of the orbits[r_1],...,[r_K]in the curve counterpart of the shape space with respect to the full distanced_F. In other words, if_1,...,_Kare the control vectors ofr_1,...,r_K, then the control vector_mean ∈ofr_meanis one representative, up to rotations, of the Fréchet mean of[_1],...,[_K]. Also, the notion of alignment proposed in <cit.> coincides with that involved in the full distance (see (<ref>)), if working with pre-shapes.

 §.§ B-spline curves


We find it is worthwhile to present our setting for B-splines, which are popular tools to construct interpolating curves (<cit.>). We provide explicit expressions for the basis functionsϕ_n. Spline curves (interpolating curves, more generally) motivate the introduction of a Hermitian productand a notion of shift configuration$̆, for which projecting on the subspace orthogonal to $̆ is akin to centering the spline curve. Readers can rely on this example to adapt our framework to other representations, such as Hermite-spline curves. 

B-spline curves have parameters called control points that are regularly spaced along a continuous parametert. A B-spline withMcontrol points holdsN = Mdegrees of freedom, as in

    ∀ t ∈ [0,1],    r(t) = ∑_n = 0^M-1[n] ϕ_n(t),

where the basis functionsϕ_ncan be obtained, for instance, from the cubic B-spline generator function

    β^3(t) = {[ |t|^3/2 - t^2 + 2/3    for  0 ≤ |t| ≤ 1; -.9cm (2 - |t|)^3/6    for  1 ≤ |t| ≤ 2;            -2.5cm 0          otherwise.;                     ].


Closed cubic B-splines curves are then obtained with

    for  n = 0, ∀ t ∈ [0,1],        ϕ_0(t) = β^3(Mt) + β^3(Mt-M),
    for  n = 1, ∀ t ∈ [0,1],        ϕ_1(t) = β^3(Mt - 1) + β^3(Mt-1-M),
    ∀ n = 2,...,M-2, ∀ t ∈ [0,1],        ϕ_n(t) = β^3(Mt-n), 
    for  n = M-1, ∀ t ∈ [0,1],        ϕ_M-1(t) = β^3(Mt + 1) + β^3(Mt + 1 - M).


One can check that1is a cubic B-spline and that it corresponds to the shift configuration=̆ (1,...,1).
We know that the temporal mean of the spline curve, defined after Definition <ref>, coincides with the product^̆* /||̆_^2 = ^̆* = ^* .
In the case of closed B-splines, it is simply the usual arithmetic mean of the control points

    r̅ = 1/M∑_m = 0^M-1[m],

which is subtracted towhen centering the configuration. Please note that we have more sophisticated expressions for closed and open Hermite-spline curves, that we detail in a future work.



 §.§ Riemannian structure of the shape space

We provide a concise description of the Riemannian structure of the shape space(our method does not use it in a crucial way). Using the same conventions as in Appendix <ref>, we use bold and italic to denote the complex and real version of a vector by∈^Nandz ∈^2N, respectively. The real inner product(z,w) ↦(^* )canonically associated to the Hermitian inner productis denoted as⟨·, ·⟩_. The shape space based on(^2N,⟨·, ·⟩_)has a Riemannian structure that is similar to the usual complex projective spaceℙ^N-2. The latter is the shape space based on(^2N,⟨·, ·⟩_can), defined for landmarks and the standard inner product <cit.>.

The Riemannian structure ofis inherited from that of, a sphere in finite dimension whose structure is well understood.
Being a smooth submanifold of(^2N,⟨·, ·⟩_),is endowed with the induced Riemannian metric. The tangent space toatzis

    T_z() = {v ∈^2N | ⟨ z, v ⟩_ = 0 },

which is equipped with the Riemannian inner product∀v_1,v_2 ∈T_z , ⟨v_1, v_2 ⟩_z := ⟨v_1, v_2 ⟩_.

The space spaceis the quotient ofby the group of planar rotationsU(1). The action ofU(1)on the Riemannian manifoldis smooth, free, and proper on, withU(1)being a Lie group acting by isometries, meaning that

    ∀θ∈ U(1), ⟨ v_1, v_2 ⟩_z = ⟨R_θ ⊙  v_1, R_θ ⊙  v_2 ⟩_R_θ ⊙   z.

The resulting quotientconsequently inherits a (unique) Riemannian structure on= / U(1)such thatΠ : ( [     →   ;  z  ↦ [] ] )is a Riemannian submersion <cit.>. 
This structure can be described as follows <cit.>. 
At each pointzof, the tangent spaceT_z()can be split into two subspaces orthogonal to each other with respect to⟨·, ·⟩_as

    T_z() = T_z([z]) ⊥⊕ H_z(),

where the vertical subspaceT_z([z]) = Ker(dπ_z) = (R_π/2  ⊙  z)is the space tangent atzto the orbit[z], while the orthogonal complement is the horizontal subspace
    H_z() = {v ∈^2N | ⟨ z, v ⟩_ = 0  and ⟨ z, R_π/2 ⊙  v ⟩_ = 0 }≃{∈̌^N  | ^* =̌ 0}.

The last identification uses the properties of the Hermitian inner product seen in Appendix <ref>. By definition, the tangent spaceT_[z]()is then identified to the horizontal subspaceH_z()through the isometrydπ_z : v ∈H_z() ↦[v] ∈T_[z]().
The Riemannian inner product onT_[z]()is then obtained as

    ⟨⟨ [v_1], [v_2] ⟩⟩_[z] := ⟨ v_1, v_2 ⟩_z,

which does not depend on the particular choice of the representative pre-shapez.

Geodesics and Riemannian distanceρ(Definition <ref>) can then be defined on. The Riemannian distanceρ_onis equal to the length of the shortest path joiningztow,

    ρ_(z,w) = arccos⟨ z, w ⟩_.

The corresponding shortest geodesic arc{(t)}that joinsztowis

    ∀ t ∈ [0,1],   (t) = 1/sin r(sin((1-t)r) z + sin(t r)w),   where  r = ρ_(z,w).

The shortest geodesic path between the shapes[]and[]is the quotient path{[(t)]}, whereis the geodesic inthat joinsz̃ ∈[]tow, wherez̃is chosen so that the corresponding (complex) pre-shapeis optimally rotated along.
The Riemannian distance in the shape space is then the length of(t), computed asarccos⟨z̃, w ⟩_= arccos((^* )). It also corresponds to the geodesic distance betweenz̃and the set[]. From an argument in the proof of Proposition <ref>, we know that, sinceis optimally rotated, it is also equal toarccos|^* | = arccos|^* |, thus resulting in the equality of Definition <ref>.


siamme