Chafik Samir

P.-A Absil

Anuj Srivastava

Eric Klassen
email: klassen@math.fsu.edu

A Gradient-Descent Method for Curve Fitting on Riemannian Manifolds *

Keywords: curve fitting, steepest-descent, Sobolev space, Palais metric, geodesic distance, energy minimization, splines, piecewise geodesic, smoothing, Karcher mean

Given data points p0, . . . , pN on a manifold M and time instants 0 = t0 < t1 < . . . < tN = 1, we consider the problem of finding a curve γ on M that best approximates the data points at the given instants while being as "regular" as possible. Specifically, γ is expressed as the curve that minimizes the weighted sum of a sum-of-squares term penalizing the lack of fitting to the data points and a regularity term defined, in the first case as the mean squared velocity of the curve, and in the second case as the mean squared acceleration of the curve. In both cases, the optimization task is carried out by means of a steepest-descent algorithm on a set of curves on M . The steepest-descent direction, defined in the sense of the first-order and second-order Palais metric, respectively, is shown to admit simple formulas.

Introduction

We are interested in the problem of fitting smooth curves to given sets of points on nonlinear manifolds. Let p 0 , p 1 , . . . , p N be a finite set of points on a Riemannian manifold M , and let 0 = t 0 < t 1 < ... < t N = 1 be distinct and ordered instants of time. The problem of fitting a smooth curve γ on M to the given points at the given times involves two goals of conflicting nature. The first goal is that the curve should fit the data as well as possible, as measured, e.g., by the real-valued function E d defined by:

E d (γ) = N i=0 d 2 (γ(t i ), p i ), (1)
where d denotes the distance function on the Riemannian manifold M . The second goal is that the curve should be sufficiently "regular" in a certain sense, e.g., the changes in velocity or in acceleration should be minimized. We are thus facing an optimization problem with two objective functions-a fitting function E d and the regularity function E s -whose domain is a suitable set of curves on the Riemannian manifold M . Curve fitting problems on manifolds appear in various applications. To cite but one example, let (I i ) i≤N be a temporal sequence of images of a 2D or 3D object motion, in which the object can appear and disappear at arbitrary times due to obscuration and other reasons. The task is to estimate the missing data and recover the motion of the object as well as possible. It is clear that focussing on the first goal (fitting the data) without concern for the second goal (regularity of the curve) would yield poor motion recovery, and that the result is likely to be improved if inherent regularity properties of the object motion are taken into account.

Previous work

One possible way of tackling an optimization problem with two objective functions is to turn it into a classical optimization problem where one of the objective functions becomes the objective function and the other objective function is turned into a constraint.

Let us first discuss the case where the fitting objective function E d is minimized under a regularity constraint. When M = R n , a classical regularity constraint is to restrict the curve γ to the family of polynomial functions of degree not exceeding m, (m ≤ N ). This least-squares problem cannot be straightforwardly generalized to an arbitrary Riemannian manifold M because the notion of polynomial does not carry to M in an obvious way. An exception is the case m = 1; the polynomial functions in R n are then straight lines, whose natural generalization on Riemannian manifolds are geodesics. The problem of fitting geodesics to data on Riemannian manifold M was considered in [START_REF] Machado | Fitting smooth paths on Riemannian manifolds[END_REF] for the case where M is the special orthogonal group SO(n) or the unit sphere S n .

The other case is when a regularity criterion E s is optimized under a constraint on E d , in which case it is natural to impose the interpolation constraint E d (γ) = 0. For example, when M = R n , minimizing the function E s defined by

E s,1 (γ) = 1 2 1 0 γ(t) 2 dt (2)
yields the piecewise-linear interpolant for the given data points and time instants (this follows from [Mil63, p. 70]), while minimizing E s,2 (γ) = 1 2 1 0 γ(t) 2 dt yields solutions known as cubic splines. (From now on, we will frequently omit the variable t in the integrand when it is clear from the context.) For the case where M is a nonlinear manifold, several results on interpolation can be found in the literature. Crouch and Silva Leite [START_REF] Crouch | Geometry and the dynamic interpolation problem[END_REF] (or see [START_REF] Crouch | The dynamic interpolation problem: on Riemannian manifolds, Lie groups, and symmetric spaces[END_REF]) generalized cubic splines to Riemannian manifolds, defined as curves γ that minimize, under interpolation conditions, the function

E s,2 (γ) = 1 2 1 0 D 2 γ dt 2 , D 2 γ dt 2 γ(t) dt, (3)
where D 2 γ dt 2 denotes the (Levi-Civita) second covariant derivative and •, • x stands for the Riemannian metric on M at x. (The subscript may be omitted if there is no risk of confusion.) They gave a necessary condition for optimality in the form of a fourth-order differential equation, which generalizes a result of Noakes et al. [START_REF] Noakes | Cubic splines on curved spaces[END_REF]. Splines of class C k were generalized to Riemannian manifolds by Camarinha et al. [START_REF] Camarinha | Splines of class C k on non-Euclidean spaces[END_REF]. Still in the context of interpolation on manifolds, but without a variational interpretation, we mention the literature on splines based on generalized Bézier curves, defined by a generalization to manifolds of the de Casteljau algorithm; see [START_REF] Crouch | The de Casteljau algorithm on the Lie group and spheres[END_REF][START_REF] Altafini | The de Casteljau algorithm on se(3)[END_REF][START_REF] Popiel | Bézier curves and C 2 interpolation in Riemannian manifolds[END_REF]. Recently, Jakubiak et al. [START_REF] Jakubiak | A two-step algorithm of smooth spline generation on Riemannian manifolds[END_REF] presented a geometric two-step algorithm to generate splines of an arbitrary degree of smoothness in Euclidean spaces, then extended the algorithm to matrix Lie groups and applied it to generate smooth motions of 3D objects. Another approach to interpolation on manifolds consists of mapping the data points onto the affine tangent space at a particular point of M , then computing an interpolating curve in the tangent space, and finally mapping the resulting curve back to the manifold. The mapping can be defined, e.g., by a rolling procedure, see [START_REF] Hüper | On the geometry of rolling and interpolation curves on S n , SOn, and Grassmann manifolds[END_REF][START_REF] Kume | Shape-space smoothing splines for planar landmark data[END_REF].

Another way of tackling an optimization problem with two objective functions is to optimize a weighted sum of the objective functions. In the context of curve fitting on manifolds, this approach was followed by Machado et al. [START_REF] Machado | Riemannian means as solutions of variational problems[END_REF] using the first-order smoothing term (2) and by Machado and Silva Leite [START_REF] Machado | Fitting smooth paths on Riemannian manifolds[END_REF] for the second-order smoothing term (3). Specifically, in [START_REF] Machado | Riemannian means as solutions of variational problems[END_REF], the objective function is defined to be

E 1 ≡ 1 2 N i=0 d 2 (γ(t i ), p i ) + λ 2 1 0 γ, γ dt,
over the class of all piecewise smooth curves γ : [0, 1] → M , where λ (> 0) is a smoothing parameter. Solutions to this variational problem are piecewise smooth geodesics that best fit the given data. As shown in [START_REF] Machado | Riemannian means as solutions of variational problems[END_REF], when λ goes to +∞, the optimal curve converges to a single point which, for certain classes of manifolds, is shown to be the Riemannian mean of the data points. When λ goes to zero, the optimal curve goes to a broken geodesic on M interpolating the data points.

In [START_REF] Machado | Fitting smooth paths on Riemannian manifolds[END_REF], the objective function is defined to be

E 2 ≡ 1 2 N i=0 d 2 (γ(t i ), p i ) + λ 2 1 0 D 2 γ dt 2 , D 2 γ dt 2 dt
over a certain set of admissible C 2 curves. The authors give a necessary condition of optimality that takes the form of a fourth-order differential equation involving the covariant derivative and the curvature tensor along with certain regularity conditions at the time instants t i , i = 0, . . . , N [MS06, Th. 4.4]. The optimal curves are approximating cubic splines: they are approximating because in general γ(t i ) differs from p i , and they are cubic splines because they are obtained by smoothly piecing together segments of cubic polynomials on M , in the sense of Noakes et al. [START_REF] Noakes | Cubic splines on curved spaces[END_REF]. It is also shown in [START_REF] Machado | Fitting smooth paths on Riemannian manifolds[END_REF]Prop. 4.5] that, as the smoothing parameters λ goes to +∞, the optimal curve converges to the best least-squares geodesic fit to the data points at the given instants of time. When λ goes to zero, the approximating cubic spline converges to an interpolating cubic spline [MS06, Prop. 4.6].

Our approach

In this paper, rather than trying to solve directly the fourth-order differential equation obtained in [START_REF] Machado | Fitting smooth paths on Riemannian manifolds[END_REF] (a feat that is not attempted in [START_REF] Machado | Fitting smooth paths on Riemannian manifolds[END_REF], except for M = R n ), we propose to search for an optimizer of the objective function using a steepest-descent method in an adequate set of curves on the Riemannian manifold M . In this section, we present the essence of our approach, and delay the mathematical technicalities until Section 2. We consider the problem of minimizing the objective function of [START_REF] Machado | Riemannian means as solutions of variational problems[END_REF], namely

E 2 : Γ 2 → R : γ → E 2 (γ) = E d (γ) + λE s,2 (γ) = 1 2 N i=0 d 2 (γ(t i ), p i ) + λ 2 1 0 D 2 γ dt 2 , D 2 γ dt 2 dt, (4)
where Γ 2 is an adequate set of curves on M to be defined in Section 2. The steepest-descent direction for E 2 is defined with respect to the second-order Palais metric v, w 2,γ = v(0), w(0) γ(0) + Dv dt (0), Dw dt (0)

γ(0) + 1 0 D 2 v dt 2 , D 2 w dt 2 γ(t)
dt.

(5)

As we shall see in Section 4, this choice of metric ensures that the gradient of E 2 at γ, represented by a vector field G along γ, admits a simple expression. This simple expression makes it possible to implement a steepest-descent algorithm on Γ 2 , where the next iterate is obtained from the current iterate γ using a line-search procedure along the path → γ on Γ 2 defined by

γ (t) = exp γ(t) (-G(t));
see Section 5. We use an Armijo backtracking procedure, but other stepsize selection methods would be suitable.

We also present a gradient-descent approach for the objective function of [START_REF] Machado | Riemannian means as solutions of variational problems[END_REF], namely

E 1 : Γ 1 → R : γ → E 1 (γ) = E d (γ) + λE s,1 (γ) = 1 2 N i=0 d 2 (γ(t i ), p i ) + λ 2 1 0 γ, γ dt, (6)
where Γ 1 is another adequate set of curves on M defined below. For E 1 , the steepest-descent direction is considered with respect to the first-order Palais metric

v, w 1,γ = v(0), w(0) γ(0) + 1 0 Dv dt , Dw dt γ(t) dt. ( 7
)
This choice confers a simple expression to the gradient; see Section 3.

Observe that the parameter λ makes it possible to balance between the two conflicting goals mentioned above: when λ is large, a higher emphasis is on the regularity condition relative to the fitting condition, whereas when λ is small, the fitting condition dominates.

The rest of the paper is organized as follows. Section 2 deals with the choice of the curve spaces Γ 1 and Γ 2 . An expression for the gradient of E 1 , resp. E 2 , is given in Section 3, resp. 4. The steepest-descent method is presented in Section 5. Numerical illustrations are given in Section 6 for M = R 2 and M = S 2 . Section 7 contains final remarks.

Preliminaries

In this section, we exploit results of Palais [START_REF] Palais | Morse theory on Hilbert manifolds[END_REF]§13] and Tromba [START_REF] Tromba | A general approach to Morse theory[END_REF]§6] to define Γ in such a way that the gradient of E with respect to the Palais metric is guaranteed to exist and be unique.

First-order case

We first consider the objective function E 1 defined in (6). Let I denote the unit interval [0, 1] and let H 0 (I, R n ) denote the set of square integrable functions from I to R n . The set H 0 (I, R n ) is a Hilbert space under pointwise operations and with the inner product •, • 0 defined by

v, w 0 = 1 0 v(t), w(t) dt,
where •, • is the standard inner product in R n . Let H 1 (I, R n ) denote the set of absolutely continuous maps γ : I → R n such that γ ∈ H 0 (I, R n ). Note that absolute continuity is equivalent to requiring that γ(t) exists for almost all t ∈ I, that γ is summable, and that

γ(t) = γ(0) + t 0 γ(s) ds.
Then H 1 (I, R n ) is a Hilbert space under the inner product •, • 1 defined by v, w 1 = v(0), w(0) + v, ẇ 0 .

(8)

This inner product belongs to a class of Riemannian structures proposed by Linnér [START_REF] Linnér | Symmetrized curve-straightening[END_REF]§3].

Let M be a closed C k+4 -submanifold of R n (k ≥ 1). Define H 1 (I, M ) to be the set of all γ ∈ H 1 (I, R n ) such that γ(I) ⊆ M . Then H 1 (I, M ) is a closed C k -submanifold of the Hilbert space H 1 (I, R n ). We set Γ 1 = H 1 (I, M ), (9)
which ensures that E 1 (6) is a well defined C k map, provided that, for all i, p i is in the image of the domain of injectivity of the exponential mapping at γ(t i ) (see Lazard and Tits [START_REF] Lazard | Domaines d'injectivité de l'application exponentielle[END_REF] for the case where the manifold is a Lie group).

The tangent space to H 1 (I, M ) at a curve γ ∈ H 1 (I, M ) is given by

T γ H 1 (I, M ) = {v ∈ H 1 (I, T M ) : v(t) ∈ T γ(t) M for all t ∈ I},
where T M denotes the tangent bundle of M . Moreover, H 1 (I, M ) is a complete C k -Riemannian manifold in the Riemannian structure induced on it as a closed C k -submanifold of H 1 (I, R n ). Note that the induced Riemannian structure on H 1 (I, M ) induced by (8) is the "extrinsic" structure given by v(0), w(0) + v, ẇ 0 , where v and ẇ are the derivatives in the sense of the embedding space R n . It thus differs from the "intrinsic" first-order Palais metric defined in (7). However, the extrinsic and intrinsic Riemannian structures are equivalent on bounded sets [Tro77, Prop. 6.1].

From this, it follows that, given γ ∈ H 1 (I, M ), the tangent space T γ H 1 (I, M ) endowed with the inner product ( 7) is a Hilbert space. The gradient of E 1 at γ is defined to be the unique

G ∈ T γ H 1 (I, M ) that satisfies, for all w ∈ T γ H 1 (I, M ), G, w 1,γ = DE 1 (γ)[w],
where DE 1 (γ)[w] denotes the derivative of E 1 at γ along w. The existence and uniqueness of G are guaranteed by the Riesz representation theorem. We will use the notation ∇E(γ) for the gradient of a function E at γ, or simply G when E and γ are clear from the context.

Second-order case

We now turn to the objective function E 2 defined in (4). Let H 2 (I, R n ) be the set of maps γ : I → R n with γ ∈ H 1 (I, R n ) and γ ∈ H 1 (I, R n ). Then H 2 (I, R n ) is a vector space under pointwise operations, and the map

Φ : R n × R n × H 0 (I, R n ) → H 2 (I, R n ) : (γ 0 , γ0 , h) → γ, defined by γ(0) = γ 0 , γ(0) = γ0 , γ(t) = h(t) for all t ∈ I, is an isomorphism. In H 2 (I, R n ), consider the inner product •, • defined by v, w = v(0), w(0) + v(0), ẇ(0) + 1 0 v(t), ẅ(t) dt.
Then Φ is an isometry and

H 2 (I, R n ) is a Hilbert space. Let M be a closed C k+4 -submanifold of R n (k ≥ 1). Define H 2 (I, M ) to be the set of all γ ∈ H 2 (I, R n ) such that γ(I) ⊆ M . Then, by restricting the proof of [Pal63, Th. 6.6] to H 2 , one obtains that H 2 (I, M ) is a closed C k -submanifold of the Hilbert space H 2 (I, R n ). We set Γ 2 = H 2 (I, M ), (10)
which ensures that E 2 is well defined. The tangent space to

H 2 (I, M ) at a curve γ ∈ H 2 (I, M ) is given by T γ H 2 (I, M ) = {v ∈ H 2 (I, T M ) : v(t) ∈ T γ(t) M for all t ∈ I}.
Given γ ∈ H 2 (I, M ), consider the mapping

Φ : T γ(0) M × T γ(0) M × H 0 (I, T γ(0) M ) → T γ H 2 (I, M )
that maps (v 0 , v0 , g) to the vector field v along γ defined by

v(0) = v 0 , Dv dt (0) = v0 , D 2 v dt 2 (t) = P t←0 γ g(t),
where P t←0 γ is the parallel transport along γ. Recall that the parallel transport is an isometry. The map Φ is an isomorphism of vector spaces between its domain and image, and it is an isometry with the obvious metric on the domain and the second-order Palais metric (5) on the image. Since the domain of Φ is a Hilbert space, its image is also a Hilbert space endowed with the inner product (5). Hence the Riesz representation theorem applies.

Gradient of E 1 in the first-order Palais metric

We derive an expression for the gradient of E 1 = E d + λE s,1 (6) over Γ 1 (9) endowed with the first-order Palais metric (7). The gradient evaluated at a curve γ involves the operations of parallel transport and covariant integral along γ.

Derivative of E d

We first give an expression for the derivative of the i th term in E d , namely,

f i : Γ 1 → R : γ → 1 2 d 2 (γ(t i ), p i ).
Let exp p denote the Riemannian exponential map at p ∈ M ; see, e.g., [START_REF] Boothby | An introduction to differential manifolds and Riemannian geometry[END_REF][START_REF] Perdigão | Riemannian geometry. Mathematics: Theory & Applications[END_REF]. Since M is a closed Riemannian submanifold of R n , it follows that M is complete (see [START_REF] Palais | Morse theory on Hilbert manifolds[END_REF]p. 326]), which means that exp p ξ exists for all ξ ∈ T p M . If q ∈ M is not in the cut locus of p, then there exists a unique minimizing geodesic α pq with α pq (0) = p and α pq (1) = q (see [dC92, corollary 13.2.8]), and we define exp -1 p (q) = αpq (0). Note that in this case, it also holds that p is not in the cut locus of q (see [dC92, corollary 13.2.7]), and we have exp -1 q (p) = -αpq (1). An expression for the derivative of f i is readily obtained from the following result.

Theorem 3.1 [START_REF] Karcher | Riemannian center of mass and mollifier smoothing[END_REF]. Let M be a complete Riemannian manifold, let p be a point of M and let q be a point of M that is not in the cut locus of p. Then the squared distance function to p is differentiable at q and we have, for all ξ ∈ T q M ,

1 2 Dd 2 (p, •)(q)[ξ] = ξ, -exp -1 q p .
Proof. This proof is essentially a restriction of the proof of [

since ċp (0, t) = 0, = α(t), -exp -1 α(t) p . Since 1 2 Dd 2 (p, •)(q)[ξ] = 1 2 d dt d 2 (p, α(t)) t=0
, the result follows. In view of this result, we have that the derivative of

f i at γ along w ∈ T γ (Γ 1 ) is Df i (γ)[w] = w(t i ), v i , where v i = -exp γ(ti) (p i ),
provided that γ(t i ) is not in the cut locus of p i . This is a mild condition, since the cut locus has measure zero [GHL04, lemma 3.96]. Finally, the derivative of E d is given by

DE d (γ)[w] = N i=0 w(t i ), v i .

Gradient of E d

The gradient of f i at γ, with respect to the first-order Palais metric (7), is the unique element g i of T γ Γ 1 such that, for all w ∈ T γ Γ 1 ,

g i , w 1 = Df i (γ)[w].
The next theorem gives an expression for g i .

Theorem 3.2. The gradient of the function

f i : Γ 1 → R : γ → d 2 (γ(t i ), p i ) evaluated at γ ∈ Γ 1 is the vector field g i along γ defined by g i (t) = (1 + t)ṽ i (t), 0 ≤ t ≤ t i (1 + t i )ṽ i (t), t i ≤ t ≤ 1 , where v i = -exp -1 γ(ti) (p i ) ∈ T γ(ti) M and ṽi is the parallel transport of v i along γ. Proof: See Appendix A.1.1.
Observe that g i is covariantly linear from 0 to t i , and is covariantly constant from t i to 1. In other words, the covariant derivative of g i is covariantly constant (ṽ i ) until t i , and it is 0 after that. Note also that ṽi (t i ) = v i .

Once we have the gradient for each of the terms in E d , the gradient of E d , under the first-order Palais metric, is simply their sum

G 1 = N i=0 g i .
(11)

Derivative of E s,1

The derivative and gradient of E s,1 (2) can be readily deduced, e.g., from [START_REF] Tromba | A general approach to Morse theory[END_REF]§6] or [KS06, Th. 1]. We give a full development here for convenience. Recall that

E s,1 (γ) = 1 2 1 0 γ(t), γ(t) dt.
Define a variation of γ to be a smooth function h : [0, 1] × (-, ) → M such that h(t, 0) = γ(t) for all t ∈ [0, 1]. The variational vector field corresponding to h is given by w(t) = h s (t, 0) where s denotes the second argument in h. Thinking of h as a path of curves in M , we define F (s) as the energy of the curve obtained by restricting h to [0, 1] × {s}. That is,

F (s) = 1 2 1 0 h t (t, s), h t (t, s) dt .
We now compute,

F (0) = 1 0 Dh t ds (t, 0), h t (t, 0) dt = 1 0 Dh s dt (t, 0), h t (t, 0) dt = 1 0 Dw dt (t), γ(t) dt, since h t (t, 0) is simply γ(t)
. Hence the derivative of E s,1 at γ along w is given by

DE s,1 (γ)[w] = 1 0 Dw dt (t), γ(t) dt.

Gradient of E s,1

In view of the above expression for the derivative of E s,1 , the following result directly follows from Section A.1.2.

Theorem 3.3. The vector field H 1 along γ that provides the gradient of the function E s,1 with respect to the first-order Palais metric satisfies the equation:

DH 1 dt (t) = γ(t), H 1 (0) = 0. ( 12
)
In the case M = R n , the gradient vector field is simply H 1 (t) = γ(t) -γ(0).

Gradient of E 1

Since E 1 = E d + λE s,1 , the gradient of E 1 follows directly from the gradients of E d and E s,1 computed above. We thus have that ∇E 1 = G 1 + λH 1 , with G 1 given (11) and H 1 given by ( 12).

4 Gradient of E 2 in the second-order Palais metric

Recall that

E 2 = E d + λE s,2 is defined on Γ 2 (10) by E 2 (γ) = 1 2 N i=0 d 2 (γ(t i ), p i ) + λ 2 1 0 D 2 γ dt 2 , D 2 γ dt 2 dt. ( 13
)
The purpose of this section is to obtain an expression for the gradient of E 2 with respect to the second-order Palais metric (5).

Gradient of E d

The derivative does not depend on the metric, in contrast to the gradient. Thus we have, as in Section 3,

Df i (γ)[w] = w(t i ), v i ,
where

f i denotes the function γ → 1 2 d 2 (γ(t i ), p i ) and v i = -exp -1 γ(ti) (p i ).
Theorem 4.1. The gradient of the function

f i : Γ 2 → R : γ → d 2 (p i , γ(t i ))
at γ ∈ Γ 2 with respect to the second-order Palais metric (5) is given by the vector field g i along γ defined by

g i (t) = (1 + t i t + 1 2 t i t 2 -1 6 t 3 )ṽ i (t) 0 ≤ t ≤ t i (1 + tt i + 1 2 tt 2 i -1 6 t 3 i )ṽ i (t) t i ≤ t ≤ 1,
where ṽi is the parallel transport of v i along γ.

Proof: See Appendix A.2.1. This gradient function is a cubic polynomial before t i and is a linear polynomial after t i . The total gradient is given by G 2 (t) = N i=0 g i (t). Another way of writing this summation is: for

t i-1 ≤ t ≤ t i , we get G 2 (t) = i-1 j=0 (1 + tt j + 1 2 tt 2 j - 1 6 t 3 j )ṽ j (t) + N j=i (1 + t j t + 1 2 t j t 2 - 1 6 t 3 )ṽ j (t). ( 14
)
4.2 Derivative of E s,2

Let γ(s, t) be a collection of curves indexed by s; for a fixed s we have a curve parameterized by t. For s = 0 that curve is simply called γ(t). Define w = ∂γ(s,t) ∂s | s=0 as the tangent vector at γ. Then DE s,2 (γ)[w] = d ds F (s) s=0 , where

F (s) = 1 2 1 0 D dt ( ∂γ ∂t ), D dt ∂γ ∂t ) dt.
Taking the derivative with respect to s:

d ds F (s) = 1 0 D ds ( D dt ( ∂γ ∂t )), D dt ( ∂γ ∂t ) dt = 1 0 [R(w, ∂γ ∂t )( ∂γ ∂t ) + D dt ( D ds ( ∂γ ∂t ))], D dt ( ∂γ ∂t ) dt,
where R is the Riemannian curvature tensor defined as:

R( ∂γ ∂s , ∂γ ∂t )(v) = D ds D dt (v) - D dt D ds (v).
(Note that the curvature tensor is sometimes defined with the opposite sign in the literature.) Since D ds ( ∂γ ∂t ) = D dt ( ∂γ ∂s ), the desired derivative at s = 0 becomes:

d ds F (s)| s=0 = 1 0 [R(w, γ)( γ) + D 2 dt 2 (w)], D dt ( γ) dt = 1 0 R(w, γ)( γ), D dt ( γ) dt + 1 0 D 2 dt 2 (w), D dt ( γ) dt. ( 15
)
This is the sought expression for DE s,2 (γ)[w].

Gradient of E s,2

We will analyze the two terms in (15) separately. The Riemannian curvature tensor has certain symmetries: for vector fields a, b, c, d along γ,

R(a, b)(c), d = -R(b, a)(c), d = -R(a, b)(d), c = R(c, d)(a), b ,
which allows us to rewrite the first term of (15) as

1 0 R( D 2 γ dt 2 (t), γ(t))( γ(t)), w(t) dt ≡ 1 0 A(t), w(t) dt .
Note that this equation defines a vector field A along the curve γ. We need a vector field H 2 with the property that H 2 , w 2 = A, w . In view of Appendix A.2.2, the solution is given by

H 2 (t) = Ĥ(t) -[-S(t) + t( Q(t) -S(t)) + 1 2 t 2 ( Q(t) -S(t)) + 1 6 t 3 S(t)], ( 16
)
where Ĥ is the four times covariant integral of A (so that it satisfies

D 4 Ĥ dt 4 (t) = A(t)) with initial conditions Ĥ(0) = D Ĥ dt (0) = D 2 Ĥ dt 2 (0) = D 3 Ĥ dt 3 (0) = 0,
and where Q and S are the parallel transport along γ of Q = D 2 Ĥ dt 2 (1) and of S = D 3 Ĥ dt 3 (1). We now consider the second term in (15), that is,

1 0 D 2 dt 2 (w), D dt ( γ) dt.
In view of Section A.2.3, this term can be written as H 3 , w 2 , where

H 3 satisfies D 2 H 3 dt 2 = D 2 γ dt 2 , H 3 (0) = DH 3 dt (0) = 0, (17)
that is, H 3 is two times covariant integral of D 2 γ dt 2 with initial conditions H 3 (0) = DH3 dt (0) = 0. In case M = R n , the two terms are simply H 2 (t) = 0 and H 3 (t) = γ(t) -γ(0)t -γ(0) for all t ∈ I.

Gradient of E 2

Combining the two gradient terms, we get the gradient of E 2 under the second-order Palais metric:

∇E 2 = G 2 + λ(H 2 + H 3 ),
where G 2 is given in ( 14), H 2 in (16), and H 3 in (17).

Steepest-descent algorithm on the curve spaces

Let E stand for E 1 (6), resp. E 2 (4), Γ for the set Γ 1 (9), resp. Γ 2 (10), of curves on the Riemannian manifold M , and let Γ be endowed with the first-order Palais metric (7), resp. second-order Palais metric (5). We propose the steepest-descent method for E described in Algorithm 1.

The algorithm creates a sequence of curves (γ k ) k=0,1,... ⊂ Γ with decreasing energy E(γ k ). The initialization step consists in choosing an arbitrary curve in Γ to be the starting curve γ 0 . Then, given the current iterate γ k , the algorithm computes the gradient ∇E(γ k ) and updates the curve to γ k+1 according to

γ k+1 (t) = exp γ k (t) (-ρ k ∇E(γ k )(t)), t ∈ I,
where ρk is a step size chosen using some step size selection rule (see, e.g., [START_REF] Bertsekas | Nonlinear programming[END_REF]). We have chosen a modified version of Armijo backtracking procedure by imposing strong Wolfe conditions [START_REF] Wolfe | Convergence conditions for ascent methods[END_REF]; see Algorithm 2. The algorithm is stopped when a certain pre-determined stopping criterion is satisfied. The criterion can be a threshold on the norm of ∇E(γ k ), for example. Whereas analyzing the convergence of steepest-descent type methods on finite-dimensional manifolds is relatively simple (see [START_REF] Absil | Accelerated line-search and trust-region methods[END_REF]), the convergence analysis of steepest-descent methods Algorithm 1 Gradient descent 1: Given a scalar ∈]0, 1[ and an initial curve γ 0 , arbitrary element of Γ; 2: k := 0; 3: repeat 4:

k := k + 1; 5:
Compute E(γ k ) and ∇E(γ k );

6:

Find the step size ρk using algorithm 2;

7: Set γ k (t) = exp γ k-1 (t) (-ρ k ∇E(γ k )(t)); 8: until ∇E(γ k ) ≤ 9: return γ := γ k Algorithm 2
Step size selection 1: Given scalars ρ 0 ∈]0, 1[, > 0 very small, 0 < σ 1 < σ 2 < 1, a function f and a descent direction q of f at x; 2: set k = 0 and set

β k = ρ 0 ; 3: until (β k ≤ ) or (f (x + β k q) ≤ f (x) + σ 1 β k ∇f (x), q and | ∇f (x + β k q), q | ≤ σ 2 | ∇f (x), q |) do 4: k := k + 1; 5:
β k = β k-1 ρ 0 ; 6: end 7: return ρ := β k on infinite-dimensional spaces is no trivial matter; see [START_REF] Smyrlis | Local convergence of the steepest descent method in Hilbert spaces[END_REF] and references therein. Analyzing the convergence of Algorithm 1 is the object of ongoing research. Nevertheless, it is reasonable to expect that the algorithm behaves like steepest-descent methods in finite dimension: the sequence of iterates γ k has a single limit (see [START_REF] Absil | Convergence of the iterates of descent methods for analytic cost functions[END_REF]) which, unless the initial curve is maliciously chosen, is a local minimizer of the objective function E. These expectations are corroborated by our numerical experiments; see Section 6.

6 Illustration on some specific manifolds:

M = R 2 , S 2
In this section we present some illustrations of our gradient descent approach to finding optimal curves. In the case of Euclidean spaces, it is sometimes possible to derive expressions for the optimal curves under E 1 and E 2 directly. In those situations, we can compare our numerical solutions to the analytical expressions, and characterize the performances. In the remaining cases, where the analytical solutions are not readily available, we will simply illustrate the results obtained using our procedures. Examples involving the analytical expressions will have M = R n and while the other cases will have M = S 2 .

Analytical solution of E 1 in R 2

As the first example we will consider the problem of finding the optimal curves under E 1 when M = R 2 . For simplicity, we will take λ = 1 in (6). This case is simple enough to seek an analytical expression as follows. Let N = 2 and the three data points be given by p 0 = (-A, 0), p 1 = (0, B), p 2 = (A, 0), at the time instants t 0 = 0, t 1 = 0.5, t 2 = 1, respectively; here A, and B be two real variables. Using the symmetry of the given points, we will note that q 0 = (-a, c), q 1 = (0, b), q 2 = (a, c) will be the control points of an intermediate curve given by the gradient descent method. Our goal is to find the values of variables a, b, and c such that the piecewise geodesic curve connecting q 0 , q 1 , q 2 is a minimum of E 1 . By computing E d and E s,1 manually we get

E 1 = 2(A -a) 2 + 2c 2 + (B -b) 2 + 4(a 2 + (b -c) 2
). The critical points are given by the equation ∇E 1 = 0, i. e. in terms of partial derivatives we have ∂E1 ∂a = ∂E1 ∂b = ∂E1 ∂c = 0. This system has only one solution given by a = 1 3 A, b = 3 7 B, and c = 2 7 B, and the minimum of E 1 is given by the piecewise geodesic curve connecting points q 1 = (-A/3, c), q 2 = (0, 3B/7), and q 3 = (A/3, 2B/7).

Shown in Figures 1((a) and (b)) are two optimal curves under E 1 obtained by our algorithm, for two different initial conditions. In each case, the green curve shows the initial condition and the black curve shows the final result obtained numerically. The red curves show the optimal obtained using the analytical solution. The coincidence of black and read curves shows the accuracy and the stability of our algorithm. In Figures 1((c) and(d)) we show the variation of the step length, and the variation of the cost function E 1 , respectively versus iterations, corresponding to the example shown in Figure 1(a).

In Figure (2) we present some additional results for R 2 , this time restricting only to our numerical solutions. These examples use a random set of points and different values of λ to demonstrate the strength of the algorithm. Each of these solutions are piecewise geodesics and the end points of the geodesic segments depend on the value of λ.

Analytical solution of E 2 in R n

Next we derive the optimal curves under E 2 for Euclidean spaces. It is interesting to note that in this case the cost function has n components, each corresponding to a coordinate in R n . In other words, the problem breaks down into n independent problems, each being one-dimensional. Therefore, it is sufficient to illustrate the analytical solution for the one-dimensional case.

To derive an analytical solution to the one-dimensional problem, we will first estblish a number of relations that this curve must satisfy and then use those relations to solve for the unknowns. We start with the fact that γ(t) = Ḧ3 (t) for all t. Therefore, γ takes the form:

γ(t) = H 3 (t) + rt + s ,
where r and s are two constants. Next, since γ is a critical point of E 2 , we have G 2 (t) = -H 3 (t) (assuming λ = 1) for all t which makes γ(t) = -G 2 (t)+rt+s. Enumerating the different conditions on γ, we obtain the following constraints.

1. Since H 3 (0) = 0, we have G 2 (0) = 0 which implies:

G 2 (0) = n j=1 v j = 0 ,
where v j is as defined in Section 4.1.

2. Also, since Ḣ3 (0) = 0, we have Ġ2 (0) = 0 which means:

Ġ2 (0) = n j=1 v j t j = 0 .
3. Finally, since we know that γ(t i ) = v i + p i , we get:

-G 2 (t i ) + rt i + s = v i + p i ,
which give us the relations: for i = 0, . . . , N

- N j=0 (1 + t j t i + 1 2 t j t 2 i - 1 6 t 3 i )v j - i-1 j=1 (1 + t i t j + 1 2 t i t 2 j - 1 6 t 3 j )v j + rt i + s = v i + p i . (18)
Rearranging this equation, we reach

n j=1 β j,i v j + rt i + s = p i ,
where

β j,i =    -[(1 -1 6 t 3 j ) + (t j + 1 2 t 2 j )t i ] j < i -[(2 + t j t i + 1 2 t j t 2 i -1 6 t 3 i )] j = i -[(1 + t j t i + 1 2 t j t 2 i -1 6 t 3 i )] j > i
Taking these three types of relations, we form a linear system of equations. We have N + 3 equations and N + 3 unknowns: show two examples of optimal curves under E2 obtained using our numerical approach (black curve) and the analytical solution (red curve). The panels in (c) and (d) plot the evolutions of E2 versus iterations for the cases (a) and (b), respectively. The curves in green are used as the initial curves for the optimization.

      1 1 . . . 1 0 0 t 0 t 1 . . . t N 0 0 β 0,0 β 1,0 . . . β N,0 t 0 1 . . . β 0,N β 1,N . . . β N,N t N 1                 v 0 v 1 v 2 . . . v N r s           =         0 0 p 0 p 1 . . . p N         . (19)
Figure 4: The optimal curves under E2 for different combinations of data points and λs.

After solving for the v j s, r and s, we can evaluate the optimal curve γ(t) = -G 2 (t) + rt + s. We present some examples for comparing the numerical solutions with this analytical solution for n = 2. In the panels Figure 3 (a) and (b), we present two examples with three points each and solve for the optimal curves under different λs. In each case, the green line shows the initial curve, the black line shows the optimal curve obtained numerically, and the red line shows the analytical solution. We have used a dotted pattern for the black curve since the two optimal curves match perfectly and one hides the other. As predicted by the theory, the optimal solution resembles a straight line when λ is sufficiently large, and an interpolating spline when λ is sufficiently small. The plots in panels (c) and (d) show the corresponding evolution of E 2 versus iterations.

In Figure 4, we present some additional examples of optimal curves (obtained using our numerical method) under E 2 in R 2 for a variety of data points and λs. Each panel in this figure shows the optimal γs for different (mostly small) values of λ but with the same data points. In each case the initial curve for the gradient process is given by the green curve.

Optimal Curves on the Unit Sphere

In this section we consider the case of M = S 2 where the analytical expressions for the optimal curves are not readily available, and we apply our numerical approach to find the solutions. In these experiments, we first generate N + 1 data points p 0 , p 1 , . . . , p N randomly on S 2 and associate them with different instants of time 0 = t 0 < t 1 < t 2 < . . . t n ≤ t 1 = 1. Then, we initialize our algorithm by an arbitrary continuous curve γ 0 ∈ Γ, and finally apply our gradient descent method to search for the optimal curve γ that minimizes E.

1. Case 1: In the case E = E 1 we apply our algorithm as described in Section 3 and examples are shown in Figure 5. Similar to the Euclidean case, the solutions are piecewise geodesic curves. Since geodesics on S 2 are arcs that lie on great circles, these optimal curves are piecewise arcs. The panels (a) and (b) show examples of optimal γ for N = 2 (three data points) and N = 3 (four data points) with λ values being 100 and 1, respectively. For λ = 100, the resulting optimal curve looks like a point. The remaining two panels (c) and (d) show several optimal curves, each corresponding to different λs, for the same set of data points. As earlier, the initial condition for the gradient descent is given by the green curve.

The Figure 6 shows two examples of the actual optimization process where the iterative updates for γ under the gradient of E 1 are shown. The process starts with the green curves as the initial conditions and the updates are shown in black. The final curves in each case are shown in red.

Case 2:

In the case E = E 2 , we need to obtain an expression for the tangent vector field A defined in Section 4.3, which involves the Riemannian curvature tensor R on M . To this end, we rely on the extrinsic expression of R given in Section B, which we particularize to the case of the sphere S 2 embedded in R 3 . The orthogonal projector onto the tangent space to S 2 at x ∈ S 2 is given by

P x = (I -xx T ), hence II(η x , ξ x ) = (D ηx P )ξ x = -η T x ξ x x and U ηx (v x ) = (D ηx P )v x = -x T v x η x ,
and thus, by (26),

R(ξ x , η x )ζ x = (η T x ζ x )ξ x -(ξ T x ζ x )η x .
Then, for the vector field A defined in Section 4.3, we have the expression

A = γ, γ D 2 γ dt 2 - D 2 γ dt 2 , γ γ. ( 20
)
Using this expression, we first integrate A(t) covariantly to determine the term H 2 of the gradient of E 2 , and then use the gradient descent method of Algorithm 1 to minimize E 2 . Shown in Figure 7 are some examples of our approach applied to different sets of points generated randomly on S 2 • If λ is very small, we have E E d . When E = E 1 (6), the optimal curve is the piecewise geodesic passing through the given points. When E = E 2 (4), the optimal curve is a piecewise cubic polynomial (in the sense of [START_REF] Noakes | Cubic splines on curved spaces[END_REF]) interpolating the given set of points when E = E 2 .

• If λ is very large, then E E s . When E = E 1 , the optimal curve shrinks to one point in M , precisely the Karcher mean of the given set of points p 0 , p 1 , . . . , p N . When E = E 2 , the optimal curve approaches the best least-squares geodesic fit to the given points.

Concluding remarks

We have addressed the problem of fitting a curve to data points on a Riemannian manifold M by means of a Palais-based steepest-descent algorithm applied to the weighted sum of a fitting-related and a regularity-related cost function. As a proof of concept, we have used the simple regularity cost function (2) based on the first derivative. We have also considered the more challenging case of the regularity cost function (3), whose derivative involves the Riemannian curvature tensor on M , and for which the optimal curves are generalized cubic splines. We have illustrated the proposed method on fitting problems in R 2 and S 2 . In future work, we will consider other nonlinear manifolds with applications in pattern recognition and image analysis.

An important feature of our approach is that the discretization takes place as late as possible in the implementation. The gradient of the cost function at a curve γ is a (continuous-time) vector field along γ expressed by means of the Riemannian logarithm, parallel transport, covariant differentiation, and covariant integrals. It is these operations that are approximated by discretization in the algorithm implementation. The advantage of using a continuous formulation is that tools from functional analysis become available. We are able to use the Palais metrics and, thus, simplify the gradient vector fields only because of this continuous formulation. An alternate approach would be to consider a discretization γ of γ using P points and discretize the function E accordingly to obtain a new objective function Ê : M P → R : γ → Ê(γ) that we would optimize on the finitedimensional product manifold M P using, e.g., a steepest-descent method described in [START_REF] Absil | Optimization Algorithms on Matrix Manifolds[END_REF]. The two approaches yield considerably different expressions for the gradient. In particular, in the approach on Ê, the gradient of the fitting term d 2 (γ(t i ), p i ) vanishes everywhere except at time t i , whereas with the approach proposed here the influence of the ith fitting term is spread along the whole curve in the expression of its gradient. Although we have not compared the results, one should expect a better performance with the approach where the discretization is delayed until the implementation step.

A.2.2 Zeroth-order L

Given a curve γ ∈ Γ 2 and a vector field A on M along γ, we seek a vector field H such that, for all w ∈ T γ Γ 2 , H, w 2 = 1 0 A, w γ(t) dt, that is,

H(0), w(0) + H (0), w (0) + 1 0 D 2 H(t) dt 2 , D 2 w(t) dt 2 dt = 1 0 A(t), w(t) dt.
We proceed by integrating by parts on the left twice, each time lowering the number of primes on w and raising the number of primes on H in the integrand:

H(0), w(0) + H (0), w (0) + H (1), w (1) -H (0), w (0) -H (1), w(1) + H (0), w(0)

+ 1 0 H (4) (t) • w(t) dt = 1 0 A(t) • v(t) dt.
From this equation, it follows immediately that H must satisfy the ODE H (4) = A, with boundary conditions H(0) + H (0) = 0, H (0) -H (0) = 0, H (1) = 0, and H (1) = 0. One can find such an H explicitly as follows: First, covariantly integrate A(t) four times to obtain a vector field H that satisfies H(4) (t) = A(t) with initial conditions H(0) = H (0) = H (0) = H (0) = 0. Note that H now satisfies the first two of the required boundary conditions on H, but not the third and fourth. So we adjust it as follows. Define Q = H (1) and S = H (1), extending each of these by parallel transport to covariantly constant vector fields Q and S along γ. Then define

H(t) = H(t) -[-S(t) + t( Q(t) -S(t)) + 1 2 t 2 ( Q(t) -S(t)) + 1 6 t 3 S(t)].
It is easy to verify that H now satisfies all four required boundary conditions. It is readily checked that the solution G satisfies G(0) = 0, DG dt (0) = 0, and D 2 G dt 2 (t) = D 2 γ dt 2 (t) for all t ∈ I.

B An extrinsic expression for the curvature tensor

The expression of the derivative of E s,2 obtained in Section 4.2 involves the Riemannian curvature tensor of the Riemannian manifold M . This tensor is well defined regardless of any embedding of the Riemannian manifold M in a Euclidean space. Nevertheless, when M is a Riemannian submanifold of a Riemannian manifold N , the curvature tensor admits an extrinsic expression in terms of the second fundamental form and of the Weingarten map, which turns out to be handy in certain cases. In this section, we present this extrinsic expression, then we work out in detail the particular case where N is a Euclidean space.

Let D η ξ denote the (Levi-Civita) covariant derivative of ξ along η on the embedding manifold M , and let D denote that derivative on N . The curvature tensor R is defined by

R(ξ, η)ζ = D ξ D η ζ -D η D ξ ζ -D [ξ,η] ζ, (23)
where ξ, η, ζ are tangent vector fields on M and [ξ, η] denotes the Lie bracket. (Observe the sign convention for R, which is not standard through the literature.) Given x ∈ M , let P x : T x N → T x M denote the orthogonal projector onto T x M . Let T ⊥ x M denote the normal space to M at x, and let P ⊥

x : T x N → T ⊥ x M denote the orthogonal projector onto the normal space T ⊥ x M . The shape operator (also called second fundamental form) is the object II defined as follows: for all x ∈ M and all ξ x , η x ∈ T x M , II(η x , ξ x ) := P ⊥ x D ηx ξ,

where ξ is any tangent vector field that extends ξ x . This definition can be found, e.g., in [O'N83,[START_REF] Chavel | Riemannian geometry[END_REF]. The Weingarten map is the object U defined as follows. For all x ∈ M , η x ∈ T x M ,

v x ∈ T ⊥ x M , U ηx (v x ) := -P x D ηx v, ( 25
)
where v is any normal vector field that extends v x . This definition can be found, e.g., in [START_REF] Chavel | Riemannian geometry[END_REF]p. 62]. Then the curvature tensor can be expressed as follows:

R(ξ, η)ζ = U ξ II(η, ζ) -U η II(ξ, ζ). (26)
Let us now assume that N is a Euclidean space. Then the projector field P can be viewed as a matrix-valued function, the shape operator admits the expression

II(η x , ξ x ) = (D ηx P )ξ x , (27)
where D now reduces to the classical derivative, and the Weingarten map takes the form

U ηx (v x ) = (D ηx P )v x . (28)
We refer to [START_REF] Absil | All roads lead to Newton: Feasible second-order methods for equality-constrained optimization[END_REF] for details. These formulas are particularly useful when the projector field P admits a simple expression.

Figure 1 :

 1 Figure 1: (a) and (b): The minimum of E1 in M = R 2 reached by the gradient descent method with respect to Palais metric using different starting curves for λ = 1, (c): the step length variation, and (d): the energy evolution versus iterations for the example shown in (a).

Figure 2 :

 2 Figure 2: The minimum of E1 in M = R 2 reached by the gradient descent method with respect to first-order Palais metric using different values of λ.

Figure 3 :

 3 Figure 3: The panels (a) and (b)show two examples of optimal curves under E2 obtained using our numerical approach (black curve) and the analytical solution (red curve). The panels in (c) and (d) plot the evolutions of E2 versus iterations for the cases (a) and (b), respectively. The curves in green are used as the initial curves for the optimization.

Figure 5 :

 5 Figure 5: Optimal curves under E1 for M = S 2 obtained by our gradient descent method with respect to the first-order Palais metric. (a): λ = 100, (b): λ = 1, (c) and (d): λ = 10, 1, 0.1, and 0.0001. In each case the green curve shows the initial condition.

Figure 6 :

 6 Figure 6: The evolution of curves under the gradient iterations for minimizing E1. (a): λ = 10 -5 , (b): λ = 100.

Figure 7 :

 7 Figure 7: Optimal curves under E2 for M = S 2 obtained by our gradient descent method with respect to the second-order Palais metric. (a): λ = 0.01, (b): λ = 10 -8 , (c) and (d): λ = 0.01, 0.001, 10 -4 , and 10 -6 . In each case the green curve shows the initial condition.

.

  The pictures in (a) and (b) show examples of optimal curves for three and five points with λs as indicated there. The remaining two panels show examples of optimal curves obtained for fixed data points under different λs. Curves in different colors are obtained by using different values of λ. The values of λ used in (c) and (d) are 0.01, 0.001, 10 -4 , and 10 -6 . As the value of λ increases, we can see the optimal curves straightening and shortening into single arcs.

Figure ( 7 )

 7 Figure (7) shows two examples of the iterative process by displaying the intermediate curves also. The initial curves are shown in green, the iterations are shown in black and the final curves are shown in red.Asymptotics on λ: Our numerical experiments corroborate the following theoretical results mentioned in Section 1.1:

Figure 8 :

 8 Figure 8: Evolutions of γ under the gradient of E2. The green curves are the initial conditions and the red curves are the final states.

A. 2

 2 .3 Second-order L Given a curve γ ∈ Γ 2 and a vector field A on M along γ, we are interested in finding a vector field G on M along γ such that, for all w ∈ T γ Γ 2 , G, w 2 = 1 0 D 2 w dt 2 , A(t) dt, that is, G(0), w(0) + G (0), w (0) + dt 2 , A(t) dt.

  Since the cut locus is closed [dC92, corollary 13.2.10], this expression is well defined for all t in a neighborhood of 0 and all s ∈ [0, 1]. Denote c p = d ds c p (s, t) and ċp = d dt c p (s, t). We know that d(p, α(t)) = c p (s, t) is independent of s.

	which still does not depend on s, thus
			=	0	1	D ds	ċp (s, t), c p (s, t) ds
			=	0	1	d ds	ċp (s, t), c p (s, t) ds
	since D ds c p (s, t) = 0 (geodesic property),
			= ċp (1, t), c p (1, t)
			We have successively
	1 2	d dt	d 2 (p, α(t)) =	1 2	d dt	c p (s, t), c p (s, t)
	which does not depend on s,
			=		D dt	c p (s, t), c p (s, t)
			=		D ds	ċp (s, t), c p (s, t)

[START_REF] Karcher | Riemannian center of mass and mollifier smoothing[END_REF] Th. 1.2]

. Let α be defined by α(t) = exp q (tξ). Consider the family of geodesics from p to α(t): c p (s, t) = exp p (s exp -1 p α(t)).

* This paper presents research results of the Belgian Network DYSCO (Dynamical Systems, Control, and Optimization), funded by the Interuniversity Attraction Poles Programme, initiated by the Belgian State, Science Policy Office. The scientific responsibility rests with its authors. This research was supported in part by AFOSR FA9550-06-1-0324 and ONR N00014-09-10664. ://www.

A Dictionary of Gradients in Palais Metrics

In a number of instances in this paper, we reach an expression for the derivative of an energy term and we want to deduce an expression for the gradient in terms of one of the Palais metrics (7), (5). In order to simplify this conversion, we derive a set of formulas for the expressions we come across in this paper.

We will use the following setup. For j = 1, 2, we are given a linear map L γ : T γ Γ j → R : w → L γ (w), and we seek an expression for G ∈ T γ Γ j such that G, w j,γ = L γ (w) for all w ∈ T γ Γ j .

A.1 First-Order Palais Metric

We first consider the case j = 1, i.e., •, • 1 is the first-order Palais metric (7) and Γ 1 is as in (9).

A.1.1 Pointwise L

For a fixed t i , let v i ∈ T γ(ti) M be an arbitrary tangent vector and let L γ (w) = w(t i ), v i γ(ti) . Our goal is thus to find a vector field g i along γ such that g i , w 1 = w(t i ), v i . We start with the left side:

Integrating the middle term by parts, we get:

Equating it with the right side, we get

Now comparing the different terms on both the sides, we get the following rules for establishing the vector field g i :

The solution is given by

where ṽi is the parallel transport of v i along γ. (Therefore, ṽi (t i ) = v i .) g i is covariantly linear from 0 to t i and is covariantly constant from t i to 1. In other words, the covariant derivative of g i is covariantly constant (ṽ i ) until t i , after that it is 0.

A.1.2 First-order L

Let γ ∈ Γ 1 , let A be a vector field along γ, and let L γ be defined by L γ (w) = 1 0 A(t), Dw dt dt. We thus seek a vector field G along γ such that

for all w ∈ T γ Γ 1 . From this expression it is clear that G must satisfy the initial condition G(0) = 0 and the ordinary (covariant) differential equation DG dt = A(t).

A.2 Second-Order Palais Metric

We now consider the case j = 2, i.e., •, • 2 is the second-order Palais metric (5) and Γ 2 is as in (10).

A.2.1 Pointwise L

For a fixed t i , let v i ∈ T γ(ti) M be an arbitrary tangent vector and let L γ (w) = w(t i ), v i γ(ti) . Hence our goal is to find a vector field g i along γ such that g i , w 2 = w(t i ), v i . Using the fundamental theorem of calculus, we have:

Therefore,

γ(t) dt = ṽi (0), w(0) γ(0) + tṽ i (t), Dw dt (t)

Applying the fundamental theorem one more time, we get:

γ(0)

Combining the last two results, we obtain

Setting it equal to w, g i 2 , we obtain the following properties for g i :

Solving for g i , we obtain:

(1 + t i t + 1 2 t i t 2 -1 6 t 3 )ṽ i (t) 0 ≤ t ≤ t i (1 + tt i + 1 2 tt 2 i -1 6 t 3 i )ṽ i (t) t i ≤ t ≤ 1 .

(22)


Sebastian Kurtek

Qian Xie

Chafik Samir

Michel Canis

Statistical Model for Simulation of Deformable Elastic Endometrial Tissue Shapes

Keywords: endometriosis, statistical shape model, simulation, elastic deformation

Statistical shape analysis plays a key role in various medical imaging applications. In particular, such methods provide tools for registering, deforming, comparing, averaging, and modeling anatomical shapes. In this work, we focus on the application of a recent method for statistical shape analysis of elastic parametrized surfaces to simulation of realistic endometrial tissue shapes. The clinical data used here contains ten magnetic resonance imaging (MRI) endometrial tissue surfaces, which are used to learn a generative shape model. We generate random samples from this model, and apply elastic semi-synthetic deformations to the randomly generated tissue shapes. This provides two types of simulated data: (1) MRI-type (without deformation) and (2) corresponding transvaginal utltrasound (TVUS) type endometrial tissue shapes, which undergo a deformation due to the transducer's pressure. The proposed models can be used for validation purposes of automatic, multimodal image registration techniques, which are crucial steps in diagnosing endometriosis.

Introduction

Endometriosis is a complex gynecological disease in which endometrial cells appear outside their usual locations in the uterine cavity [START_REF] Brosens | Diagnosis of endometriosis: pelvic endoscopy and imaging techniques[END_REF]. The main symptoms depend on the site of active endometriosis and are influenced by hormonal changes. This disease affects approximately 10% of women in the reproductive age group and may cause chronic pelvic pain, severe dysmenorrhea, infertility, rectal bleeding and digestive problems. Endometriosis can be found in the pelvic cavity region, specifically in the pelvic peritoneum and pelvic organs. Currently there is no etiologic cure for endometriosis, but it can be treated in a variety of ways including pain medication, hormonal treatments, and laparoscopic surgery in severe cases. An accurate diagnosis must be made to obtain important medical information. Diagnosis and surgery planning in endometriosis are often improved by analyzing shapes of organs and tissues. Recent advances in medical imaging offer increasingly detailed information on typical anatomical structures. However, there is a lack of validation techniques for automatic image registration strategies, especially for multimodal images. In many medical applications, real data can only be extracted manually by an expert, and then used to validate image processing algorithms. Indeed, scarcity of data for evaluation results in restricted studies. In this paper, we present a new statistical framework to generate realistic data that can be used as ground truth when dealing with deformability of endometrial cells. Standard methods to assess an accurate diagnosis use multiple modalities including trans-vaginal ultrasound (TVUS) and magnetic resonance imaging (MRI). However, some limitations due to non-localized endometrial lesions or its infiltration in other organs cannot be directly avoided. An interesting solution is to statistically analyze shapes of real clinical data and provide enough random or simulated samples to validate the TVUS to MRI registration step; registration of these two modalities is key for fusing complementary information for diagnostic purposes [START_REF] Yavariabdi | Contour-based TVUS-MR image registration for mapping small endometrial implants[END_REF][START_REF] Samir | Elastic shape analysis of cylindrical surfaces for 3D/2D registration in endometrial tissue characterization[END_REF].

There are many approaches to generate synthetic data for validation of medical image processing methods [START_REF] Commowick | Incorporating priors on expert performance parameters for segmentation validation[END_REF][START_REF] Sotiras | Deformable medical image registration: A survey[END_REF]. First, one could use physical image phantoms with known shapes. However, phantoms are limited to specific information with restricted variability and are very costly if adapted to complex medical cases. Second, virtual organs can be simulated using controlled numerical models. These approaches are usually based on parametric models with good approximations and computational cost. But, they are limited to controlled deformations governed by the model. Thus, ignoring the variability of anatomical structures limits the realism of these simulations. Finally, one can utilize tools from shape analysis to characterize large amounts of natural variability. In this paper, we adopt the square-root normal field (SRNF) representation [START_REF] Jermyn | Elastic shape matching of parameterized surfaces using square root normal fields[END_REF] of cylindrical surfaces to provide a set of tools, including parallel transport, exponential map, and geodesics, needed for statistical analysis of endometrial tissue shapes.

Clinical Data and Contributions:

This study was carried out using real data from ten patients who have small endometrial implants in the pelvic area. For each patient, MRI was used to examine their pelvic organs. First, MRI slices which include both the endometrial implant and its neighboring organs were selected for each patient. MR images used for these experiments had an average size of 400 ⇥ 400 ⇥ 5mm 3 with a voxel resolution of 0.5 ⇥ 0.5 ⇥ 5mm 3 . Second, the soft tissue organs (i.e. bladder, uterus, ovary, rectum) and the implant were segmented by an expert. Finally, the endometrial tissue data was represented using a cylindrical surface parameterization, which was constructed from a set of 2D MR contours. Additionally, corresponding 2D TVUS images were segmented to provide deformed endometrial tissue curves. An example of the reconstruction process is shown in Figure 1. All of the surfaces in our dataset, reconstructed using a basis defined in the following section, are shown in Figure 2 (the original data is shown in [START_REF] Kurtek | Statistical shape model of elastic endometrial tissue surfaces[END_REF]). There is a lot of variation in this data, and thus, parsimonious shape models are very important in this application.

Some preliminary results of this study were presented in a recent conference paper [START_REF] Kurtek | Statistical shape model of elastic endometrial tissue surfaces[END_REF]. The novel contributions of the present article are (1) random sampling from a Gaussian model using the exponential map instead of a linear approximation, and (2) simulation of semi-synthetically deformed endometrial tissue shapes. The first contribution is theoretical in nature, although it can also have practical implications. In particular, we extend the tools proposed in [START_REF] Xie | Parallel transport of deformations in shape space of elastic surfaces[END_REF] to apply to cylindrical surfaces and use shooting geodesics to develop random sampling from a Gaussian endometrial tissue shape model. The second contribution allows for simulation of TVUS-type, deformed endometrial tissue shapes via parallel transport on the shape space of cylindrical surfaces. The rest of this paper is organized as follows. In Section 2, we define a mathematical framework for shape analysis of cylindrical surfaces. Section 3 describes tools for statistical modeling of endometrial tissue shapes and a procedure for generating semi-synthetic endometrial tissue deformations. Finally in Section 4, we close with a brief summary.

Mathematical Framework

Let F be the space of all smooth embeddings of a cylinder in R 3 , where each such embedding defines a parametrized surface f : S 1 ⇥ [0, 1] ! R 3 . Let be the set of all boundary-preserving diffeomorphisms of S 1 ⇥ [0, 1]. For an endometrial tissue surface f 2 F, f represents its re-parameterization. As shown in previous work, it is inappropriate to use the L 2 metric for statistical shape analysis of parameterized surfaces, because does not act on F by isometries [START_REF] Jermyn | Elastic shape matching of parameterized surfaces using square root normal fields[END_REF][START_REF] Kurtek | Parameterization-invariant shape statistics and probabilistic classification of anatomical surfaces[END_REF]. Thus, we utilize the square-root normal field (SRNF) representation of surfaces and the corresponding Riemannian metric proposed in [START_REF] Jermyn | Elastic shape matching of parameterized surfaces using square root normal fields[END_REF]; additional tools for statistical analysis were given in [START_REF] Xie | Parallel transport of deformations in shape space of elastic surfaces[END_REF]. However, most of the previous work was limited to spherical surfaces. Thus, in the following we provide a summary of these methods and particularize them to apply to cylindrical endometrial tissue surfaces.

Let s = (u, v) 2 S 1
⇥ [0, 1] be a coordinate system on the unit cylinder. Then, n(s) = @f @u (s) ⇥ @f @v (s) denotes a normal vector to the surface f at the point f (s). The SRNF representation of surfaces is defined using a mapping Q :

F ! L 2 as Q(f )(s) = n(s) |n(s)| 1/2 .
The space of all SRNFs is a subset of L 2 and it is endowed with the natural L 2 metric. The differential of the map Q, denoted by Q ⇤,f , is used to define the corresponding Riemannian metric on F as g [START_REF] Kurtek | Statistical shape model of elastic endometrial tissue surfaces[END_REF] for a detailed expression). It is now easy to verify that the re-parameterization group acts on F by isometries under the pullback metric g F ; this property is necessary to define a parameterizationinvariant shape analysis framework. Furthermore, this metric is automatically invariant to translation because it is based on partial derivatives of f , w 1 and w 2 only. Scaling variability can be removed by rescaling all surfaces to have unit area. With a slight abuse of notation, we redefine F as the space of all unit area surfaces, which forms the pre-shape space in our analysis.

F (w 1 , w 2 ) = g L 2 (Q ⇤,f (w 1 ), Q ⇤,f (w 2 )), w 1 , w 2 2 T f (F) (see
In order to define tools for computing parallel transport and geodesics, we must first define an orthonormal basis on F, B = {b 1 , b 2 , . . . }. Any surface can be expressed with respect to this basis as f (s) =

P 1 k=1 ↵ k b k (s). Then f = (↵ 1 , ↵ 2 , . . . ) 2 R 1
forms an alternative representation of f . We define the orthonormal basis B in two steps. First, we define two basis sets on S 1 and [0, 1] separately as

B S 1 = {sin(n 1 u), 1 cos(n 1 u), 1|n 1 = 1, . . . , N 1 , u 2 [0, 2⇡]}, and B [0,1] = {sin(2⇡n 2 v), 1 cos(2⇡n 2 v), 1 v, v|n 2 = 1, . . . , N 2 , v 2 [0, 1]}.
Then, to form the basis set B we take all products of the elements in B [0,1] and B S 1 , and orthonormalize using Gram-Schmidt under the L 2 metric. As a result, we obtain a total of N orthonormal basis elements; we have truncated this basis to a finite number despite F being an infinite dimensional space. Figure 2 shows the reconstruction of the original data using this basis; the surfaces are shown in the same order as in [START_REF] Kurtek | Statistical shape model of elastic endometrial tissue surfaces[END_REF].

We are interested in modeling shapes of endometrial tissue surfaces. Thus, we must also remove the variability in our data due to rotation and re-parameterization. This is also known as the registration process, where one-to-one correspondence is determined across two (or multiple) endometrial tissue surfaces. Registration can be performed in the SRNF space as follows. For a rotation matrix O 2 SO(3) and a surface f 2 F, the SRNF of the rotated surface Of is Oq. Also, for a reparameterization 2 , the SRNF of a re-parametrized surface f is (q, ) = (q ) q J , where J is the determinant of the Jacobian of . In order to regis- ter endometrial tissue surfaces with respect to rotation and re-parameterization, we first define an equivalence class of an SRNF as [q] = {(Oq, )|O 2 SO(3), 2

L(F ⇤ ) = 0.5334
}. Each equivalence class represents an endometrial tissue shape uniquely. Then, the registration problem can be stated as

(O ⇤ , ⇤ ) = arginf (O, )2SO(3)⇥ kq 1 (Oq 2 , )k 2 .
This optimization problem is solved iteratively. First, one fixes and searches for an optimal rotation over SO(3) using Procrustes analysis. Then, given this rotation, one searches for an optimal re-parameterization over using a gradient descent algorithm [START_REF] Jermyn | Elastic shape matching of parameterized surfaces using square root normal fields[END_REF]. In the following, we let

f ⇤ 2 = O ⇤ (f 2

⇤

).

To compute geodesics between shapes of endometrial tissue surfaces, we utilize a tool called parallel transport. We first give a brief summary of what parallel transport accomplishes and refer the reader to [START_REF] Xie | Parallel transport of deformations in shape space of elastic surfaces[END_REF] for more details. Recall that we have defined an orthonormal basis, B, for the vector space F. Suppose we are interested in transporting a deformation vector field v 2 T F (0) (F) along a path F (t) connecting two optimally registered surfaces

F (0) = f 1 and F (1) = f ⇤ 2 .
Since we can express F (t) and v as F (t) ⇡

P N k=1 ↵ k (t)b k and v ⇡ P N k=1 a k (0)b k , we want to define the vector field Y (t) ⇡ P N k=1 a k (t)b k along F (t) such that Y (0) = v and DY dt = 0.
In other words, we want to obtain a vector field Y (t) that has the same "direction" and magnitude as the initial vector field Y (0). The coefficients of this vector satisfy an appropriate differential equation with the initial condition

P N k=1 a k (0)b k ⇡ Y (0)
. Next, we describe two procedures for computing geodesic paths and distances between endometrial tissue surfaces.

Shooting Method (Exponential Map): Given a surface f and a deformation v 0 , we can evaluate the exponential map exp f (tv 0 ) = F ⇤ (t), t = [0, 1] using m discrete segments. First, we initialize

F ⇤ (0) = f and v(0) = v 0 . For the ⌧ th segment, given v( ⌧ 1 m ) and F ⇤ ( ⌧ 1 m ), we first set F ⇤ ( ⌧ m ) = F ⇤ ( ⌧ 1 m ) + 1 m v( ⌧ 1 m ). Then, we compute the parallel transport of v( ⌧ 1 m ) from F ⇤ ( ⌧ 1 m ) to F ⇤ ( ⌧ m
) and denote it by v( ⌧ m ). We repeat this procedure for ⌧ = 2, . . . , m. The geodesic distance in this case is given by the length of the initial deformation vector v 0 :

L(F ⇤ ) = q g F (v 0 , v 0 ).
Path Straightening: The energy of a path F on the shape space of surfaces is given by E

(F ) = R 1 0 g F ( d dt F (t), d dt F (t))dt.
It is a standard result that the critical points of E are geodesics. The gradient of E can be computed using parallel transport as follows. First, compute the vector field u(t) defined as the covariant integral of the gradient rE = w, the path F can be updated as F (t) = F (t) ✏rE, where ✏ > 0 and small. This is done until convergence resulting in the geodesic path F ⇤ . The geodesic distance is given by the length of this path:

L(F ⇤ ) = R 1 0 q g F ( d dt F ⇤ (t), d dt F ⇤ (t))dt.
In Figure 3, we show an example of a geodesic between endometrial tissue shapes computed using path-straightening.

Statistical Shape Model of Endometrial Tissue Shapes

We begin by defining an intrinsic mean shape under the proposed metric, called the Karcher mean. Let {f 1 , f 2 , . . . , f n } 2 F denote a sample of endometrial tissue surfaces. Also, let F ⇤ i denote a geodesic path between a surface f and a surface f ⇤ i that was optimally registered to f . This geodesic is computed using path-straightening. Then, the sample Karcher mean is given by f = argmin f 2F

P n i=1 L(F ⇤ i ) 2 , where F ⇤ i (0) = f and F ⇤ i (1) = f ⇤ i .
A gradient-based approach for finding the Karcher mean is given in [START_REF] Dryden | Statistical Shape Analysis[END_REF].

Once the sample Karcher mean has been computed, the evaluation of the Karcher covariance, which captures the observed variability around the mean, is performed as follows. First, we optimally register all surfaces in the sample to f , resulting in {f ⇤ 1 , . . . , f ⇤ n }. Next, we find the shooting vectors {⌫ 1 , . . . , ⌫ n } from f to each surface in {f ⇤ 1 , . . . , f ⇤ n }. We perform principal component analysis (PCA), using Gram-Schmidt under g F , to generate an orthonormal basis {B j |j = 1, . . . , k}, k  n, of the observed {⌫ i }. We project each ⌫ i onto B using ⌫ i ⇡ P k j=1 c i,j B j , where c i,j = g F (⌫ i , B j ). Now, each original surface can be represented using given data. If u 2 R k corresponds to a principal singular vector of K, then the corresponding shooting vector is given by P k j=1 u j B j ; this vector is mapped to a surface f using the exponential map. In Figure 4 (top), we display the path traced by following the two main directions of variation from 1 standard deviation (blue) to +1 standard deviation (red) around the Karcher mean (green). First, note that the computed mean is a nice representative of the real data. The observed endometrial tissue surfaces have local convex and concave structures and the Karcher mean is of similar structure. Second, these paths provide a natural set of representative deformations present in our data and an efficient summary of observed variability. Figure 4 (bottom) provides a different visualization of the first three principal directions of variation. The left panel provides the deformation vector fields on the mean surface corresponding to each positive principal direction. The right panel is a visualization of the pointwise magnitude of this vector field. This allows for easy identification of areas undergoing highest amounts of deformation.

c i = {c i,j }. Then, K = 1 n 1 P n i=1 c i c T i 2 R
We validate our model using random sampling. For this purpose, we utilize the Gaussian distribution defined in the tangent space at the mean endometrial tissue surface. A random tangent vector can be generated using v = the corresponding principal singular vector and B j is a basis element. One can then obtain a sample from the Gaussian distribution using the exponential map

f rand = exp f (v).
To simulate semi-synthetic deformations of the simulated elastic endometrial tissues, we first define a diffusion model with ten degrees of deformation on the Karcher mean; each deformation describes density dynamics in soft tissue undergoing diffusion. We parallel transport the semi-synthetic deformation vectors from the Karcher mean to each of the random samples generated by our statistical shape model. We then shoot a geodesic path in that direction to obtain the deformed random sample. As a consequence, the surfaces can be deformed in very complicated manners with varying magnitudes. A random point on the Karcher mean surface is selected as the center of the deformation, since the TVUS probe can move freely in any direction. Then, the deformation vector field on the rest of the surface is obtained by obeying the same boundary and initial conditions to share common spatial and temporal evolutions; these comprehensive deformations are learned to simulate tissue deformation and probe orientation. In Figure 5(a) we display ten examples of randomly sampled endometrial tissue shapes from the proposed Gaussian model. We also show a deformed version of these random samples in (b) based on the semi-synthetic deformations described above and displayed in panel (c) as a deformed cylinder. The generated random samples and their deformed versions naturally resemble the given data, which is a desirable property of our model. This visual evaluation was performed on 99 random samples.

Model Evaluation: Specificity, generality and compactness are three common evaluation criteria of a shape model. Compactness refers to the amount of variance in the shape model. While we do not compare the proposed model to any other methods in current literature, we point out that our model is guaranteed to be compact due to its invariance to re-parameterization. This idea was shown for curves and surfaces in [START_REF] Kurtek | Parameterization-invariant shape statistics and probabilistic classification of anatomical surfaces[END_REF][START_REF] Kurtek | Statistical analysis of manual segmentations of structures in medical images[END_REF]. Specificity refers to the ability of a shape model to represent only valid shapes, while generality quantifies the ability of the model to describe unseen shapes. Because there is a clear trade-off between these two measures, one would ideally like to strike a balance between them. In order to evaluate the proposed endometrial tissue shape models we utilize the statistical framework for analyzing shapes of curves proposed in [START_REF] Srivastava | Shape analysis of elastic curves in Euclidean spaces[END_REF]. We use curves rather than surfaces for the evaluation because the only real data available for the TVUS deformed endometrial tissues are manually segmented 2D curves. Thus, we first extract a random level curve in proximity to the center of the applied deformation from each simulated, deformed endometrial tissue shape. For consistency, we evaluate the non-deformed shapes in the same manner. In this case, we also extract a random level curve from the given MRI data. After these steps, we have 99 MRI-type and TVUS-type curves coming from the simulated endometrial tissues and 10 MRI-type and TVUS-type curves coming from real data. We then use the framework in [START_REF] Srivastava | Shape analysis of elastic curves in Euclidean spaces[END_REF] to perform the evaluation separately for the MRI and TVUS shape models. First, we estimate the mean and covariance structure for the simulated data. We perform PCA and retain the two principal directions of variation. Then, we project the simulated data and the real data onto these two directions, resulting in a two-dimensional plot shown in Figure 6 (left). Note that we could have also estimated the principal directions of variation using the real data, but due to the very small number of observations we decided to use the simulated data instead. It is clear from this plot that both shape models do a fairly good job at exploring the entire shape space of endometrial tissues satisfying generality. It is also clear that we do not stray very far from the real data, which is related to specificity. As a second step, we compute the pairwise distance matrix after combining the simulated and real data. Then, we perform multidimensional scaling (MDS) to visualize the similarity in the entire dataset. We display the result as a two-dimensional plot in Figure 6 (right). As in the PCA results, we see evidence that our shape model is both specific and general.

We also sought evaluation from a radiologist with expertise in gynecology. Here, the radiologist only evaluated the non-deformed (MRI) endometrial tissue shape model. We first presented the radiologist with geodesic paths from the estimated Karcher mean to each of the ten given observations; this provided an idea of shape deformations generated via our model. Next, the radiologist examined 50 randomly sampled shapes from the proposed Gaussian model. He was then asked to provide ten scores (one for each original observation; 1=bad, 2=acceptable, 3=good, 4=very good), based on the plausibility of the geodesic deformation, and how well the original data is represented in the generated random samples. The three main criteria considered by the radiologist, which are closely related to the generality and specificity of the model, were the consistency of the global shape based on his general experience and on the given data, the deformability of the tissues, and the deformation/distance from the Karcher mean. The scores provided by the radiologist were one 2, three 3s, and six 4s (mean=3.5, variance=0.5), which support our claim that the proposed generative model produces valid endometrial tissue shape samples.

Summary

We have defined a statistical shape model of deformable, elastic endometrial tissue shapes. The proposed model efficiently captures variability in the given clinical data and produces natural random samples resembling MRI endometrial tissues.

The described methodology allows for registration, deformation and efficient statistical modeling of endometrial tissue shapes. We utilize an additional elastic deformation model and parallel transport to generate deformed endometrial tissues as would be observed in TVUS images. The proposed framework can be used in simulation studies where a large number of observations is needed to validate automatic, multimodal image registration techniques. In our future work, we plan on utilizing this model for such validation.

Fig. 1 .

 1 Fig. 1. An example of constructing a cylindrical surface using MR image slices.

Fig. 2 .

 2 Fig. 2. Endometrial tissue surface data from [7] reconstructed using the basis B.

Fig. 3 .

 3 Fig. 3. Example of a geodesic path (and distance) between endometrial tissue surfaces.

ddtFFig. 4 .

 4 Fig. 4. Top: Two main directions of variation in the given data displayed as a path from 1 standard deviation (blue) to +1 standard deviation (red) around the Karcher mean (green). Bottom: Principal direction deformation vector fields (left) and their point-wise magnitudes (right) for the three principal directions of variation in the given data.

Fig. 5 .

 5 Fig. 5. (a) Randomly sampled shape from the Gaussian model resembling MRI data. (b) Random sample after additional deformation resembling TVUS data. (c) Deformation applied to the random sample displayed on a perfect cylinder.

SFig. 6 .

 6 Fig. 6. Principal component analysis (left) and multidimensional scaling (right) based evaluation of generality and specificity of the proposed shape model for non-deformed (MRI) and deformed (TVUS) endometrial tissue shapes. (Simulated data=blue; Real data=red)