The framework presented here, described generally in the following Section, makes a distinction between the growing body, a mathematical object, and the material it contains at each growth stage, which is a subset of another mathematical object-the material manifold .

We use a global geometric framework to continuum mechanics (see [START_REF] Epstein | Differentiable Manifolds and the Principle of Virtual Work in Continuum Mechanics[END_REF][START_REF] Segev | Forces and the Existence of Stresses in Invariant Continuum Mechanics[END_REF]) in which the collection of configurations of a system is given the structure of a differentiable manifold, whose tangent bundle contains the generalized velocities and whose cotangent bundle contains the generalized forces. Here, the collection of various material contents of the growing body (various material structures)-the content space-is an infinite dimensional manifold. Thus, the generalized forces associated with the growth processes are elements of the cotangent bundle of the content space. Similarly, the configuration space of the growing body contains all the configurations in space of all the material contents that the body may possess.

The basic mathematical objects, those describing the growing body, material manifold, a material content, and a configuration of the growing body in space, are defined in Section 2. Section 3 considers the content space and the configuration space of the growing body. In particular, the configuration space of the growing body is assumed to have the structure of a fiber bundle. A desirable feature in a theory of growing bodies is the possibility to decompose a generalized velocity of the growing body into two components: a rate of change of the material content and a material velocity field. Such a decomposition will induce a similar decomposition of forces into components associated with growth and components associated with mechanical power. This decomposition, whose geometrical counterpart is a connection on the configuration space (when regarded as a fiber bundle), is considered in Section 4.

Section 5, based on [START_REF] Segev | On Smoothly Growing Bodies and the Eshelby Tensor[END_REF], models the growth of organisms, i.e., growing bodies that have identifiable elements. This theory is analogous to the theory of volumetric growth (see [START_REF] Skalak | Analytical Description of Growth[END_REF]6]). Section 6, based on [START_REF] Segev | On Symmetrically Growing Bodies[END_REF], presents a theory capable of modeling the growth of bodies due to phase transition. This theory is obtained from the theory of organisms by the requirement of invariance under the action of the group of diffeomorphisms of the growing body. Section 7 presents an example where the material manifold is a not a Euclidean space. The situation described in Section 7 can be interpreted as the growth of wood in the cross section of a tree trunk. Section 8 formalizes the surface growth of Skalak et al. [START_REF] Skalak | Analytical Description of Growth[END_REF] used to model growth in bones, horns, spiral shells and trees. Section 9 generalized the situation of Section 8 to allow variable age of the growing bodies.

Stress theory is not discussed here; however, the relevant results of [START_REF] Segev | Forces and the Existence of Stresses in Invariant Continuum Mechanics[END_REF][START_REF] Segev | On the Consistency Conditions for Force Systems[END_REF][START_REF] Segev | Kinematics of and forces on nonmaterial interfaces[END_REF][START_REF] Segev | Force Theory for Multiphase Bodies[END_REF] concerning the representation of forces by stresses may be used to construct stress theories for the various examples considered here.

The constructions of the various infinite dimensional manifold structures use results of Kijowski & Komorowski [START_REF] Kijowski | A Differentiable Structure in the Set of All Bundle Sections over Compact Subsets[END_REF], Komorowski [START_REF] Komorowski | A Geometrical Formulation of the General Free Boundary Problems in the Calculus of Variations and the Theorems of E. Noether Connected with Them[END_REF], Binz & Fischer [START_REF] Binz | On the Manifold of Embeddings of a Closed Manifold[END_REF], Michor [START_REF] Michor | Manifolds of Differentiable Mappings[END_REF] and Hamilton [START_REF] Hamilton | The Inverse Function Theorem of Nash and Moser[END_REF].

GROWING BODIES

In this Section we describe the basic mathematical objects used in the formulation of the theories of growing bodies that we present in the following.

The notion of a growing body is used in order to extend continuum mechanics that deals traditionally with bodies of a fixed material structure so it can be used in situations in which the material structure of the body varies. This changing of material structure reflects the growth of the body.

While traditional continuum mechanics considers the body and the physical space as basic objects, theories of growing bodies replace the "body" by two objects: the growing body Band the material manifold M . We will use the term simple body when we refer to a regular body B of continuum mechanics.

It is assumed that during the process of growth, in all the various growth stages, the growing body retains certain identifiable properties. For example, one may assume that the topology of the growing body is retained during a process of growth. The mathematical object that models the growing body is chosen in such a way that it reflects those properties that remain invariant and identifiable throughout the process of growth. Thus for example, if for a certain theory the topology is the only property that remains invariant throughout processes of growth, then B will be a topological space. Henceforth, it is assumed that the growing body has a differentiable structure.

The material manifold M is conceived as the collection of virtual material points. This notion is somewhat analogous to the material universe (see Noll [19] and Truesdell [START_REF] Truesdell | A First Course in Rational Continuum Mechanics[END_REF]). Thus, a point X in M may be a part of the body in a certain growth stage or it may not. As suggested by the terminology, it is assumed that M is a smooth manifold.

A growth stage of the body is described mathematically by a mapping c : B → M. The image of the mapping c is the collection of material points that the body contains at the growth stage described by c and accordingly, we will refer to c as a content of the growing body. Thus, the image of a content is a simple body. The class of mappings that are admitted as contents should reflect the invariant properties of the growing body. For the example of a body whose topology is conserved in all its growth stages, a content will be a continuous mapping. In the following we will be more restrictive than this example and it will is assumed henceforth that the images of contents and simple bodies are smooth compact manifolds with boundaries.

A simple interpretation of these notions may be offered if one identifies the material manifold with the space manifold and the images of contents of the growing body are identified with reference configurations of simple bodies in space. This interpretation can be carried one step further if we identify the growing body B with a simple body B and consider the images of the contents of B as various reference configurations of the simple body B.

Since for each content c of the growing body, B = c{B} is a simple body, we may consider configurations of B in the physical space. For simplicity, we consider here only the case where the physical space is modeled by 3 , and in accordance with most treatments of continuum mechanics, we will model configurations of simple bodies in space by embeddings. Thus, we define a configuration of a growing body κ to be a configuration in space of the image of some content c of the growing body.

With the interpretation of a content of the growing body as a reference configuration of a simple body, a configuration of the growing body may be interpreted as a deformation superimposed on the reference configuration.

It is customary in continuum mechanics to consider the action of groups on simple bodies and on the space manifold. The action of a group on a body is usually related to symmetry of the material properties of the body and the action of a group on the space manifoldusually, the action of the Euclidean group on 3 -is generates balance laws. In the case of theories of growing bodies, the actions of groups on the growing body and material manifold should be considered.

In the following, we will present an example for the action of a group G B on the growing body manifold. We will denote this action by Ψ B so that

Ψ B : G B × B → B
is a smooth mapping. Similarly to the case of a group action on a simple body, the action Ψ B describes invariance of the "properties" of the growing body points under class of certain transformations. In addition, one may consider the action

Ψ M : G M × M → M
of a group G M on the material manifold.

Theories of growing bodies may be classified according to the invariance requirements that are satisfied. That is under what group action Ψ B on the growing body and under what group action Ψ M on the material manifold is the theory invariant. For example, if one assumes that M has a metric structure, it is natural to require that contents that differ by an isometry of M be identified. (With the interpretation of a content as a reference configuration, this requirement means identifying reference configurations of the body that differ by rigid body motions.) Such a requirement of invariance will lead to balance of content forces.

THE CONTENT SPACE AND THE CONFIGURATION SPACE OF A GROWING BODY

We will study the various theories from a global geometric point of view and so we make the following definitions. The content space C is the collection of all contents of the growing body in the material manifold, i.e., C = {c | c : B → M is a content}. It is assumed that the content space is a differentiable manifold. Clearly, as a manifold of mappings, C is infinite dimensional. When considering various theories of growing bodies, one has to show for each instance that C possesses indeed a structure of a manifold.

Given a content c, one may consider the configuration space Q c{B} of its image, the simple body B = c{B}, in the physical space. That is, 3 ), where, "Emb" denotes the collection of embeddings of the first argument in the second.

Q c{B} = Emb(B,
Thus, the configuration space of the growing body is

Q B = c∈C Q c{B} .
We note that Q B has a natural projection mapping π C :

Q B → C that assigns c to every κ in Q c{B} .
The following results will be of use when considering the structure of the various configuration spaces. The collection Emb(X , Y) is an open subset of the Frechet manifold C ∞ (X , Y) of smooth mappings of X in Y, for a compact X (see [START_REF] Michor | Manifolds of Differentiable Mappings[END_REF][START_REF] Hamilton | The Inverse Function Theorem of Nash and Moser[END_REF]). Thus, the tangent space T i Emb(X , Y) to Emb(X , Y) at the embedding i can be identified with the tangent space to

C ∞ (X , Y) at i, i.e., with {u ∈ C ∞ (X , T Y) | τ Y • u = i},
where τ denotes the tangent bundle projection of the corresponding manifold. Since, 3 is a trivial manifold, it follows that for any configuration κ of c{B} in space, T κ Q c{B} is isomorphic with C ∞ (c{B}, 3 ). We conclude that for every c in C, the fiber π -1 C (c) is a manifold. In [START_REF] Epstein | Differentiable Manifolds and the Principle of Virtual Work in Continuum Mechanics[END_REF][START_REF] Segev | Forces and the Existence of Stresses in Invariant Continuum Mechanics[END_REF][START_REF] Segev | On the Consistency Conditions for Force Systems[END_REF] the kinematics and force theory of simple bodies was presented from a global geometrical point of view. The configuration space of a simple body in space was presented and generalized velocities were defined as elements of the tangent bundle. The results quoted in the previous paragraph indicate that as expected, such generalized velocities are the virtual velocity fields of continuum mechanics. Generalized forces, defined as elements of the cotangent bundle, were shown to generalize forces on simple bodies in continuum mechanics. In particular, under a natural choice of topology for the configuration space, forces where shown to be represented by stresses. This discussion explains the terminology material velocity field for an element v of T κ Q c{B} and a simple body force for an element f m in T * κ Q c{B} , where c is a given content.

We will consider only theories of growing bodies for which the configuration space possesses a differentiable structure. For such theories it is possible to apply the framework of [START_REF] Epstein | Differentiable Manifolds and the Principle of Virtual Work in Continuum Mechanics[END_REF][START_REF] Segev | Forces and the Existence of Stresses in Invariant Continuum Mechanics[END_REF][START_REF] Segev | On the Consistency Conditions for Force Systems[END_REF] so that generalized velocities are defined as elements of the tangent bundle of the configuration space and generalized forces are defined as elements of the cotangent bundle. Thus, one may formulate kinematics and force theory for growing bodies.

Similarly, since C was assumed to have a smooth structure, generalized velocities and forces pertaining to the growth of the body may be defined as elements of the tangent and cotangent bundles of C respectively. Thus, we will refer to an element ċ ∈ T c C as a content rate and to an element f c ∈ T * c C as a content force. Clearly, the content rate and content force are velocity and force only in the generalized sense of the word. The use of such forces is motivated by the works of Eshelby, [START_REF] Eshelby | Energy Relations and the Energy-Momentum Tensor in Continuum Mechanics[END_REF], and Gurtin (e.g., [START_REF] Gurtin | Multiphase Thermomechanics with Interfacial Structure 3. Evolving Phase Boundaries in the Presence of Bulk Deformation[END_REF][START_REF] Gurtin | The Characterization of Configurational Forces[END_REF][START_REF] Gurtin | The Nature of Configurational Forces, Archive for Rational Mechanics and Analysis[END_REF]) who uses the terms accretive forces and configurational forces in his works on multiphase bodies.

The mapping π C , our assumptions concerning the availability of differential structure for both Q B and C, and the earlier discussion concerning the structure of π -1 C (c), c ∈ C, will be supplemented by the requirement that Q B has the structure of a fiber bundle. The availability of a connection on the fiber bundle will be related to our ability to decompose generalized velocities and forces pertaining to the growing body into the corresponding variables pertaining to growth and simple bodies. We will describe the problem in some detail in the next section.

Note that with the foregoing framework, the interpretation of contents as reference configurations of a simple body may be formulated mathematically as follows. We assume that we have a given section s : C → Q B of the bundle structure of Q B . Such a section associates with every content c of the growing body a configuration s(c)-the "reference configuration"-of the growing body in space.

CONNECTION-INDUCED DECOMPOSITIONS

A desirable feature that one would like a theory of growing bodies to have is the possibility to decompose a generalized velocity into a growth rate and a velocity of a simple body. Similarly, it is desirable that forces be decomposable into content forces and simple body forces. Such decompositions are made possible by a connection on the bundle π C : Q B → C. This Section describes some of the details pertaining to such decompositions.

The mapping τ Q B : T Q B → Q B assigns the configuration κ of the growing body to every generalized velocity of the growing body at that configuration. In addition, we note that the tangent mapping

T π C : T Q B → T C
assigns to every generalized velocity the corresponding content rate. Moreover, we can identify generalized velocity fields that are actually material velocity fields by the condition that the content rate associated with them vanishes, i.e., generalized velocity fields κ that satisfy the condition T π C ( κ) = 0. Geometrically, such purely material generalized velocity fields are termed vertical vectors, the subspace V κ Q B ⊂ T κ Q B containing the vertical vectors is usually called the vertical subspace and the collection of all vertical subspaces is the vertical subbundle V Q B of T Q B . Nevertheless, we do not have an invariant mapping that assigns material velocity fields to generalized velocities.

A connection on Q B is specified by means of a family of linear injections

Γ κ : T c C → T κ Q B , κ ∈ Q B , c = π C (κ) satisfying, T π C •Γ κ = 1
, where 1 is the identity mapping. The mapping Γ κ will be referred to as the connection mapping. The connection mapping induces a mapping

∆ κ : T κ Q B → V κ Q B for every κ ∈ Q B by ∆ κ = 1 -Γ κ • T π C .
We will refer to ∆ κ as the vertical projection and, indeed, one can verify that the image of ∆ κ consists of the vertical vectors. Similarly, we will refer to Image{Γ κ } = Kernel(∆ κ ) as the horizontal distribution. Using ι to denote the inclusion mapping, the situation is illustrated in the following commutative diagram.

V κ Q B ι ---→ T κ Q B Tκπ C ---→ T c C V κ Q B ∆κ ← ---T κ Q B Γκ ← ---T c C
Thus, (T π C , ∆ κ ) provide the required decomposition of generalized velocities.

Once a connection is given, one can take duals of the foregoing mappings and decompose forces by (Γ κ * , ι * ). The situation is described by the following commutative diagram.

(V κ Q B ) * ι * ← ---(T κ Q B ) * (Tκπ C ) * ←----(T c C) * (V κ Q B ) * (∆κ) * ---→ (T κ Q B ) * (Γκ) * ---→ (T c C) *

ORGANISMS

The theory of organisms describes growing bodies whose elements are identifiable throughout the growth. It is natural to assume that for a living organism we are able to identify various parts of the organism irrespective of the growth. In accordance with continuum mechanics we require that this ability does not end at a certain scale or size. This requirement is reflected mathematically in the following definitions. Definition 5.1. An organism B is a three dimensional compact manifold with boundary that can be embedded in 3 . Definition 5.2. The material manifold of an organism is a three dimensional Euclidean space M with tangent space V. Definition 5.3. A content c of an organism is an embedding of the organism in the material manifold.

Clearly, the requirement that contents are embeddings means that not only the growing body points are identifiable but in addition there is a differentiable structure that remains invariant throughout the growth.

Proposition 5.1. The configuration space of an organism has the structure of a trivializable fiber bundle.

Proof. Given a configuration κ in π -1 C (c) for c ∈ C, set, e : B → 3 by e = κ • c. We will refer to e as the extent of the growing body at the configuration κ. Consider the mapping c,e). This natural mapping defines a global fiber bundle chart on Q B .

Φ : Q B → C × Emb(B, 3 ) given by Φ(κ) = (π C (κ), κ • π C (κ)) = (
Under this global chart a generalized velocity is represented by (c, e, ċ, ė) where ċ ∈ T c C is a vector field on B valued in the tangent space V to the material manifold M , and ė ∈ T e Emb(B, 3 ) is a vector field on B valued in 3 .

The natural global fiber bundle chart Φ induces immediately a decomposition of generalized velocities for which the connection mapping Γ κ is given by T c Φ(Γ κ ( ċ)) = ( ċ, 0). Here, T c Φ denotes the tangent of Φ at c and it follows that ∆ κ (T c Φ -1 ( ċ, ė)) = T c Φ -1 (0, ė). In other words, the representative of the vertical projection ∆ κ acts as the projection on the second component in the product C × Emb(B, 3 ). From the definition of Φ it follows that ∆ κ ( κ) = ė • c -1 . With this connection a generalized velocity κ belongs to the horizontal distribution if ė = 0 so that the growing body points remain momentarily stationary in space. This means that if the body grows but at the same time it is deformed so that the identifiable parts remain in the same location in space the corresponding generalized velocity will be horizontal.

A different connection can be defined on Q B such that the generalized velocity corresponding to the situation describe above will no longer be horizontal and the vertical projection will have a clear kinematical meaning.

Proposition 5.2. The mapping Γ κ given by

Φ(Γ κ ( ċ)) = Dκ • c( ċ)
defines a connection on Q B for which the vertical component of a generalized velocity is the material velocity field

v(X) = ė • c -1 (X) -Dκ(X)( ċ • c -1 (X)), X ∈ c{B}.
Proof. One can verify that the last equation holds by differentiating the relation κ(t)(X) = e(t) • c -1 (t)(X), X ∈ c{B}, with respect to the t parameter. In addition, using the definition of ∆ κ in terms of Γ κ one can show that the expression for v is indeed the result of the vertical projection.

If f B : T κ Q B → is a force on the growing body, the global chart Φ allows a natural representation of f B in terms of a pair of functionals

(f c , f e ) ∈ C ∞ (B, V) * × C ∞ (B, 3 ) * , in the form f B ( κ) = f c ( ċ) + f e ( ė),
where ċ and ė are the representatives of the generalized velocity κ. It is natural to refer to f c as a content force and to f e as an extent force.

Corresponding to the decomposition of generalized velocities into vertical and horizontal components, one has a decomposition of generalized forces. The decomposition induced by the first connection mentioned is virtually identical to the decomposition using the natural chart except for the fact that the domain of definition of the vertical component is c{B} and not B. We consider therefore the decomposition induced by the material velocity connection. Thus, a force f B acting on the growing body at the configuration κ may be represented by

(f a , f m ) ∈ C ∞ (B, V) * × C ∞ (c{B}, 3 ) * , in the form f B ( κ) = f a ( ċ) + f m (v),
where, v is the material velocity field corresponding to κ. Naturally, we will refer to f m -a simple body force-as the material component of the force. Note that although f a belongs to the same space as f c they are different. The power expanded by f B on κ will be equal to the power expanded by f c on ċ if ė vanishes. Thus, intuitively we can say that f c performs work on both the content rate and the rate in which the material deforms.

The relation between the components of forces are given by the following equations.

f a = f c + f e • (Dκ • c) f m = f e • c * ,
where, c * is the pullback mapping of the vector fields defined on c{B} into vector fields defined on B, so for example,

c * (v)(ξ) = v(c(ξ)), v : c{B} → 3 , c * (v) : B → 3 , ξ ∈ B.

INVARIANCE UNDER THE GROUP OF DIFFEOMORPHISMS

In this section we use the same setting as in the case of organisms but now we consider a theory that is invariant under the action of the group of diffeomorphisms of B, i.e., we consider the case where G B is the group of diffeomorphisms of B, Diff(B). The mechanical situation that such a theory describes is completely different than the theory of organisms. Here, the invariance under the group of diffeomorphisms implies that the growing body points are no longer identifiable. Thus, the only property that remains invariant throughout growth is the differential topological structure of the growing body. As such, this theory is aimed at describing phenomena such as phase transition for the case where a single phase is considered. Definition 6.1. We say that two embeddings c 1 , c 2 have the same

shape if there is a diffeomorphism d ∈ Diff(B) such that c 2 = c 1 • d.
Clearly, the condition that two embeddings have the same shape defines an equivalence relation on Emb(B, M ). The equivalence classes are the shapes of the growing body and two embeddings have the same shape if and only if they have the same image. Thus, the image of a shape is well defined and in the sequel we will often identify a shape with its image. This motivates the term "shape" and the definition above. We conclude that if we wish to construct a theory of growing bodies that is invariant under the action of the group of diffeomophisms on B we have to admit shapes of B in M as contents rather then embeddings as in the previous Section. We will use χ to denote a typical shape and we will use S to denote the collection of shapes. Definition 6.2. The content space for the Diff(B)-invariant growing body theory is S = Emb(B, M )/ . The configuration space for the Diff(B)-invariant growing body theory is

Q B = χ∈S Emb(Image(χ), 3 ).
In order that this theory satisfies the postulates of Section 3 we have to show that the collection of shapes has a differentiable structure. This result, given first in [START_REF] Kijowski | A Differentiable Structure in the Set of All Bundle Sections over Compact Subsets[END_REF] can be also found in [START_REF] Binz | On the Manifold of Embeddings of a Closed Manifold[END_REF][START_REF] Michor | Manifolds of Differentiable Mappings[END_REF][START_REF] Hamilton | The Inverse Function Theorem of Nash and Moser[END_REF]. Proposition 6.1. The collection of shapes S has a structure of a Frechet manifold. The tangent space T χ S can be identified with C ∞ (∂χ), the collection of smooth functions defined on the boundary of χ.

Proof. We will only outline the construction of charts on the space of shapes. The rest of the proof is available in the references cited above. Let n be the unit normal field to ∂χ and for u in C ∞ (∂χ), consider the mapping s u : ∂χ → M given by s u (X) = X + u(X)n.

It may be shown that for a sufficiently small u, the image of s u is a two dimensional submanifold of M that is diffeomorphic to ∂χ. Thus, the image of s u is the boundary of a shape χ . In addition, for two distinct u 1 , u 2 , Image(s u 1 ) = Image(s u 2 ). In other words, once χ is given with its normal field, there is a neighborhood U of χ such that the mapping

φ : U → C ∞ (∂χ), so that φ(χ ) is the unique function in C ∞ (∂χ) with Image(s u ) = χ ,
is a chart on S. We will refer to this chart as the chart centered at χ. Remark 6.1. We note that the chart in S constructed above and the identification of T χ S with C ∞ (∂χ) depended on our choice of the normal vector field rather than any other vector field that is transversal to ∂χ. While the choice of the normal is natural here, in general, any other vector field may be used to construct a chart. Moreover, it is not necessary to use straight lines in the construction. Any family of parametrized curves, transversal to ∂χ-a tubular neighborhood -may be used. Such constructions may be needed if the material manifold M does not have a metric structure. In such a general case, the basic principle of the construction still applies but the the tangent space T χ S is identified with the space of smooth sections C ∞ T M/T (∂χ) , where T M/T (∂χ) is the quotient vector bundle-normal bundle.

The next result that we quote pertains to the manifold structure of Q B . Proposition 6.2. The configuration space for the Diff(B)-invariant growing body theory has a structure of a fiber bundle 3 ) at any shape χ can be identified with Emb(B, 3 )-the space of extents.

π C : Q B → C, whose fiber π -1 C {χ} = Emb(χ,
Proof. We will outline the construction of a fiber bundle chart in a neighborhood of a content χ ∈ S. We construct a diffeomorphism between Emb(χ , 3 ) and Emb(χ, 3 ) for χ in a neighborhood of χ. Such a diffeomorphism can be constructed by artificially deforming χ to the shape χ . Mathematically this is obtained by generating a diffeomorphism δ χ : χ → χ that depends smoothly on χ . Such a mapping δ will be referred to as a dragging of the domain χ, and its existence in a neighborhood of χ can be proved. Thus, π -1

C {χ } = Emb(χ , 3 ) can be identified with

π -1 C {χ} = Emb(χ, 3 ) by κ → κ = κ • δ χ , κ : χ → 3 , κ : χ → 3 .
Remark 6.2. The manifold Emb(B, 3 ) is diffeomorphic to Emb(χ, 3 ). However, unlike the theory of organisms, there is no natural identification of π -1 C {χ} with Emb(B, 3 ) as χ is not a particular embedding but an equivalence class.

The tangent space T κ Q B can be identified using a fiber bundle chart with 3 ). This way, a generalized velocity κ is represented by a real function on the boundary of χ and an extent rate. The first component χ is interpreted as the normal velocity field of growth. We recall that the first component is the representative of T π C ( κ) and it depends on our choice of the unit normal field, or in the general case, choice of a particular tubular neighborhood. The second component, the extent rate, depends in addition on our choice of dragging of the domain χ and the particular shape χ of B.

C ∞ (∂χ) × C ∞ (B,
The next step in the construction of a Diff(B)-invariant growing body theory is the consideration of a connection on Q B . Proposition 6.3. Using the unit normal field, there is a connection on Q B such that the vertical component of a generalized velocity κ is the material velocity field v.

We omit the proof of the Proposition. For the details, see Segev, Fried & deBotton [START_REF] Segev | Force Theory for Multiphase Bodies[END_REF] and Hamilton [START_REF] Hamilton | The Inverse Function Theorem of Nash and Moser[END_REF]. If X ∈ Interior(χ) and B(t) is the motion in Q B whose tangent at t = 0 is the generalized velocity κ, then

v(X) = d dt κ(t)(X) t=0 makes sense because X ∈ π C • κ(t) for t in a neighborhood of t = 0.
The values of v on the boundary of χ are obtained as limits. An additional connection is also suggested by Hamilton [START_REF] Hamilton | The Inverse Function Theorem of Nash and Moser[END_REF] whose mechanical application we do not consider here. Forces on a Diff(B)-invariant growing body can be represented using a fiber bundle chart by means of a pair (f S , f e ), where f S ∈ C ∞ (∂χ) * is interpreted as a force normal to the boundary of χ and f e ∈ C ∞ (B, 3 ) an extent force. Thus,

f B ( κ) = f S ( χ) + f e ( ė),
where ( χ, ė) are the representatives of κ.

Using the decomposition induced by the connection, we have a representation of forces in the form

f B ( κ) = f S ( χ) + f m (v), f S ∈ C ∞ (∂χ), f m ∈ C ∞ (χ, 3 ) * ,
where, χ = T π C ( κ), and v is the material velocity field.

A NON-EUCLIDEAN MATERIAL MANIFOLD

In this Section we present an example in which the material manifold M is not a Euclidean space. For the sake of illustration we consider a two dimensional example but the extension to three dimensions is straightforward.

Consider the case where B is a two dimensional annulus. We denote by ∂B 1 and ∂B 2 the inner and outer components of the boundary of the annulus, respectively. We set M to be the two dimensional cylinder S 1 × L, where S 1 denotes the circle and L denotes a one dimensional oriented Euclidean space-a line-together with a specific orientation. We will treat M as a one dimensional vector bundle π : M → S 1 over the circle, and much of the following immediately generalizes to the case of a general vector bundle.

Note that an embedding of the annulus into the cylinder can be in one of two forms. For the first form, the boundary of the annulus can be contracted to a point on the cylinder and in the second form it is impossible to contract the image of the boundary to a point, and ∂B 1 , ∂B 2 are the images of two sections s 1 and s 2 of π, respectively. Note also that for an embedding of the second kind, there are two possible orientations of the image of B in M . For one orientation, s 1 (X) > s 2 (X) for each X on S 1 and for the second orientation s 1 (X) < s 2 (X). Simply put, we define a content to be the image of a mapping that wraps the annulus on the cylinder like a sleeve so that the outer boundary of the annulus remains always on one particular side of the image. Definition 7.1. A content of the growing body is a shape χ of embeddings B → M such that the images of the boundary components ∂B 1 and ∂B 2 correspond to sections s 1 and s 2 of π, respectively, with s 2 (X) > s 1 (X), for all X on S 1 . Remark 7.1. Note that we can also set M = S 1 × L 2 with the bundle structure π : S 1 × L 2 → S 1 , induced by the projection on the first factor, so that a content can be identified with a section of π with s 2 > s 1 . Definition 7.2. A configuration of the growing body B at the content χ is an embedding κ : χ → 2 that preserves the orientation of the annulus.

The next Proposition follows immediately from the definition of a content and the fact that the pullback of the vector bundle π : M → S 1 onto the images of the two components of the boundary serve as tubular neighborhoods.

Proposition 7.1. The content space 2 . Given a content χ by s 1 and s 2 any shape χ = (s 1 , s 2 ) in a neighborhood of χ in C may be represented by (s 1 -s 1 , s 2 -s 2 ) which is in a neighborhood of the zero element of C ∞ (S 1 ) 2 . Here, we identify the tangent space T X L with so the difference of two smooth sections can be identified with an element of C ∞ (S 1 ).

C = {(s 1 , s 2 )|s j ∈ C ∞ (π), j = 1, 2, s 2 (X) > s 1 (X), for all X ∈ S 1 }, where, C ∞ (π) is the Frechet space of smooth sections of π, is an open subset of C ∞ (π)
Clearly, the method of dragging of the domain can be used in order to provide Q B with a differentiable structure. A configuration of the growing body κ will be represented under a fiber bundle chart on Q B by a pair (χ, e), where, χ = (s 1 , s 2 ) ∈ C ∞ (S 1 ) 2 , and e : B → 2 is an orientation preserving embedding. A generalized velocity κ is represented under a chart by the pair ( χ, ė), χ = ( ṡ1

, ṡ2 ) ∈ C ∞ (S 1 ) 2 , ė ∈ C ∞ (B, 2
). Thus, a force f B is represented under a chart by a pair (f S , f e ), where f S = (f shps 1 , f S 2 , f S j 2 being real distributions over the circle, and

f e ∈ C ∞ (B, 2 ) * .
The fact that M is not the Euclidean plane does not affect the construction of the connection outlined in the previous Section using the material velocity field as the vertical component of the generalized velocity. Using the connection, generalized velocities decompose by

κ → ( χ, v), χ ∈ C ∞ (S 1 ) 2 , v ∈ C ∞ (χ, 2 ), generalized forces decompose by f B → (f S , f m ), f S = (f S 1 , f S 2 ), f S j ∈ C ∞ (S 1 ) * , f m ∈ C ∞ (χ, 2 ) * , and finally f B ( κ) = f S 1 ( χ1 ) + f S 2 ( χ2 ) + f m (v).
Remark 7.2. In this example, motivated by the notion of surface growth of Skalak et al. [START_REF] Skalak | Analytical Description of Growth[END_REF], using a cylinder rather than 2 to model the material manifold allows unbounded growth on both components of the boundary of the annulus, even at the inner boundary. If one uses 2 or 2 -{(0, 0)} to model M , the metric properties of the material manifold will not be represented adequately.

Remark 7.3. The addition of material at the inner boundary of the annulus is meant to model the growth of the bark of a tree trunk at the cambium, while there is no generation of material at s 2 . In order to model the generation of wood at the inner side of the cambium one should consider contents that have two components c = (χ 1 , χ 2 ). Here, χ 1 describing the outer growth is defined as above and χ 2 is a shape of hte circle bounded by ∂B 1 in 2 -an + valued function on ∂B 1 . Alternatively, one should add an additional dimension to the cylinder representing M . For example, we can use the fiber bundle

π : S 1 × L 2 × + → S 1
to model the material manifold. Shapes of the annulus in S 1 × L 2 × + that are images of sections of the type (s 1 , s 2 , s 3 ), s 2 > s 1 , of this fiber bundle will model contents. The value of the third component represents the growth of wood-the circle bounded by ∂B 1 . A description of surface growth that is closer to the treatment of [START_REF] Skalak | Analytical Description of Growth[END_REF] is given in the following Section.

SURFACE GROWTH

Skalak et al. [START_REF] Skalak | Analytical Description of Growth[END_REF] use the notion of surface growth to model growth patterns in bones (particularly the skull), horns, spiral shells, and trees. The characteristics of surface growth in the work of Skalak et al. are as follows. The material points in the image of a content are parametrized by points on a surface and a parameter τ . For a particular X in the image of a content, the parameter τ represents the time elapsed since X was added to the body. The parameter τ varies between zero and a maximum value t, the age of the growing body. Growth may occur on both sides of the surface, possibly in two different directions.

We now turn a formal geometrical description of surface growth. We first consider growth on one side of the surface only and extend the situation to growth on two sides growth later in the Section. Definition 8.1. The growing body for the theory of surface growth is

B = Σ × [0, 1],
where Σ is a compact two dimensional submanifold without boundary of 3 to which we will refer as the growth surface. The material manifold M is a three dimensional Euclidean space with tangent space V. A content of the growing body is an isotopy of the growth surface in M , i.e., an embedding c : B = Σ × [0, 1] → M. Remark 8.1. By the definition, the material points in the images of contents have identifiable surface and "time" parameters, and as such, it is an example of an organism with a somewhat different interpretation of the growing body points. The fact that we consider the fixed nondimensional "time" interval [0, 1] means that the various contents of the growing body are as growth stages of the growing body at one fixed "age". Remark 8.2. Consider a two dimensional example (M is 2-dimensional) in which Σ = S 1 . We interpret Σ as a cross section of the cambium of a tree and for a content c in C, we interpret Image(c) as the cross section of a tree trunk having a fixed nondimensional age one. Growth of wood on one side of the cambium is considered only, i.e., the growth of wood rather than the bark. We note that with this model, the trunk has a cavity because of the fact that Image(c) is an annulus. The rings of the trunk are the surfaces c{Σ × {τ }}, τ ∈ [0, 1], and so, the points on the cross section are parametrized by the "time" they were created. The content space is interpreted as the collection of sections of tree trunks of a fixed age corresponding possibly to various kinds of trees and various growth conditions. We note that in the previous example, the growth of the trunk in the various radial dimensions need not be identical because of the invariant parametrization by the points on Σ. If we require that the theory is invariant under the action of the diffeomorphisms or rotations on S 1 , we obtain growth patterns that are symmetrical.

The fact that we consider a special case of the theory of organisms implies that all the properties of that theory hold here, in particular, the differentiable structures and the connection.

In order to model growth on two sides of the growing surface, the previous framework should be modified as follows.

Definition 8.2. The growing body for the theory of two sided surface growth is

B = B 1 ∪ B 2 , B 1 = Σ × [-1, 0], B 2 = Σ × [0, 1],
where Σ is a growth surface. A content of the growing body in M is a continuous injection c : B → M such that the restrictions

c 1 = c| B 1 : B 1 → M, c 2 = c| B 2 : B 2 → M, are smooth embeddings.
This definition makes contents for the theory of two sided surface growth identical in structure to composite body configurations discussed in [START_REF] Segev | Force Theory for Multiphase Bodies[END_REF]. The mathematical structure of B here is that of a composite body with the two phases B j , and interface Σ × {0}. Roughly, contents are "embeddings" of B in M whose derivative suffers a jump across the interface. The configuration space of a composite body is considered in [START_REF] Segev | Force Theory for Multiphase Bodies[END_REF] and the results apply to two sided surface growth as follows.

Proposition 8.1. The content space C for the theory of two sided surface growth is an open subset of

F M = {w ∈ C 0 (B, M ) | w j = w |B j ∈ C ∞ (B j , M ), j = 1, 2}.
Thus, C is a differentiable manifold with tangent space

F V = {u ∈ C 0 (B, V) | u j = u |B j ∈ C ∞ (B j , V), j = 1, 2}.
We note that the jump in the derivative of contents at τ = 0 implies that the differential structure of B is different than that of B = c{B}. Thus, while a configuration of the growing body in space is by definition a smooth embedding κ : c{B} → 3 , for some c ∈ C, the corresponding extent e = κ • c : B → 3 is no longer smooth but has a jump in the derivative at τ = 0, in contrast with the situation for the theory of organisms used for one sided surface growth. The structure of the configuration space should be modified accordingly.

Proposition 8.2. The configuration space for a two sided surface growth has a differentiable structure of a trivializable fiber bundle, such that

Φ : Q B → C × F 3 , given by Φ(κ) = (π C (κ), κ • π C (κ)) = (c, e),
where,

F 3 = {u ∈ C 0 (B, 3 ) | u j = u |B j ∈ C ∞ (B j , 3 ), j = 1, 2}, is a global natural trivialization.
One would also like to use the material velocity field in order to construct a connection on the configuration space as in the case of the theory of organisms. However, recalling the expression v(X) = ė • c -1 (X) -Dκ(X)( ċ • c -1 (X)), X ∈ c{B}, for the material velocity field, we note that because of the discontinuity in the derivative Dκ, the velocity field suffers a jump. Hence, the velocity field is not a member of the vertical tangent subspace

T κ (π -1 C (c)) = C ∞ (c{B}, 3 )
and the construction of a connection cannot be followed. This situation is analogous to the theory of multiphase bodies (see [START_REF] Segev | Force Theory for Multiphase Bodies[END_REF]).

SURFACE GROWTH WITH A VARIABLE AGE

The mathematical framework suggested in the previous Section was intended to model surface growth in the case where all bodies had a fixed age-that corresponding to the dimensionless time parameter τ = 1. While a material point was parametrized by the τ -parameter specifying the time it was created, the maximum value of the time parameter-the age of the body-was fixed for all bodies. In this Section we wish to extend the framework suggested for surface growth so that the contents of growing bodies may have different identifiable ages. For the sake of simplicity, we restrict the discussion to one sided growth. Definition 9.1. For a growth surface Σ, the t-old generation is

C t = Emb(Σ × [0, t], M ), t ∈ + .
The content space for a variable age surface growth theory is

C = t∈ + C t .
The mapping π t : C → + that assigns t to every element of C t will be referred to as the age projection.

The construction here is somewhat similar to the situation with the configuration space over the manifold of shapes in a Diff(B)-invariant theory and similar tools should be used for the construction of the bundle structure.

Proposition 9.1. The content space for a variable age surface growth theory has the structure of a trivializable fiber bundle π t : C → + .

Proof. To construct a global fiber bundle chart we construct a diffeomorphism of C t with C 1 for t ∈ + . The dragging of the domain Once the content space has been defined, the structure of the configuration space Q B = c∈C Emb(Image(c), 3 ) should be considered. Note that we can write

Q B = t∈ + c∈Ct
Emb Image(c), 3 .

Thus, setting

Q B t = c∈Ct
Emb Image(c), 3 to be the collection of configuration in space of the t-old generation, the configuration space can be written as the union of configuration spaces of fixed age surface growth theories,

Q B = t∈ + Q B t .
Proposition 9.2. The configuration space for a variable age surface growth theory has two fiber bundle structures. It is a trivializable fiber bundle π C : Q B → C, and a trivializable fiber bundle

π t • π C : Q B → +
where the fiber over t in + is Q B t .

Proof. We first consider fiber bundle charts for π t •π C . Let t ∈ + , then, the fiber bundle structure of π t implies that there is a diffeomorphism C t to C 1 given by c → c • δ -1 t . Hence, Next consider c in C t . We construct a diffeomorphisms of Q Image(c) with Emb(Σ × [0, 1], 3 ). Here, we only have to generalize the extent e = κ • c used for organisms (or surface growth with a fixed age), into

Q B t =
e = κ • c • δ -1 t ∈ Emb(Σ × [0, 1], 3 ), c ∈ C t , κ ∈ Q Image(c) .
Thus, surface growth with variable age provides an example of an iterated fiber bundle structure for the configuration space.

  δ t : Σ × [0, t] → Σ × [0, 1],takes the very simple form(ξ, τ ) → (ξ, τ t ), so that c → c = c • δ -1 t , c ∈ C t , c ∈ C 1 , is the required diffeomorphism.

3 is

 3 diffeomorphic to Q B 1 because for c ∈ C t , c ∈ C 1 , c = c • δ -1 t we have Image(c ) = Image(c).

Acknowledgements. This work was partially supported by the Paul Ivanier Center for Robotics Research and Production Management, Ben-Gurion University.


Pierre-Francois Marteau

Sylvie Gibet

Clément Reverdy

Down-Sampling coupled to Elastic Kernel Machines for Efficient Recognition of Isolated Gestures

In the field of gestural action recognition, many studies have focused on dimensionality reduction along the spatial axis, to reduce both the variability of gestural sequences expressed in the reduced space, and the computational complexity of their processing. It is noticeable that very few of these methods have explicitly addressed the dimensionality reduction along the time axis. This is however a major issue with regard to the use of elastic distances characterized by a quadratic complexity. To partially fill this apparent gap, we present in this paper an approach based on temporal down-sampling associated to elastic kernel machine learning. We experimentally show, on two data sets that are widely referenced in the domain of human gesture recognition, and very different in terms of quality of motion capture, that it is possible to significantly reduce the number of skeleton frames while maintaining a good recognition rate. The method proves to give satisfactory results at a level currently reached by state-ofthe-art methods on these data sets. The computational complexity reduction makes this approach eligible for real-time applications.

I. INTRODUCTION

During the past decade, gesture recognition has been a very active research field that has evolved in terms of improving motion capture technology and recognition methods, mostly based on machine learning techniques. Recently, the availability of low-cost consumer technology, often associated with game consoles, has helped to democratize the use of motion sensors, not only in the context of interactive video game, but also in various frameworks using gestural interaction. Hence high quality databases, built from expensive motion capture (mocap) devices requiring specific expertise, exist today alongside lower quality databases, i.e. containing noisier and less accurate data provided by new inexpensive sensors that require very little expertise. Therefore, heterogeneous databases of captured motion of various qualities are available to the scientific community, and comparing the robustness and generalization of recognition algorithms on these diverse motion databases is highly challenging. Beyond the quality of the recognition, the complexity of the algorithms and their response time is also a major issue, especially for real-time interaction.

In the context of the recognition of isolated gestures from motion captured data, we present in this paper a novel method that improves the performance of classical support vector machines when used with regularized elastic kernels. We also show how the temporal dimensionality reduction, associated with such elastic kernels significantly improves the efficiency of the algorithm and the recognition scores.

II. MOTION CAPTURED DATA AND SEQUENCE OF SKELETAL POSES

We focus on human motion data captured by various camera-based sensors (infrared marker-tracking system with high resolution, or webcam-style system). The data is uniformly preprocessed so that the captured data is finally reconstructed as a set of 3D-trajectories of the skeleton joints determined from the positions of markers captured on a real actor underlying more or less accurately the skeleton of the actor who produced the movement. The identification of the skeleton model from captured data is achieved through a mappingoptimization process such as the ones described in [START_REF] De Aguiar | Automatic learning of articulated skeletons from 3d marker trajectories[END_REF], [START_REF] O'brien | Automatic joint parameter estimation from magnetic motion capture data[END_REF], or [START_REF] Shotton | Real-time human pose recognition in parts from single depth images[END_REF]. The techniques based on a skeleton model hence convert 3D sensor data into Cartesian or angular coordinates that define with various accuracies the state of the joints over time.

Figure 1 presents two skeletons reconstructed from two very distinct capture systems. On the right, the skeleton is reconstructed from data acquired via the Microsoft Kinect, on the left from the Vicon-MX device used by the Max Planck Institute to produce the HDM05 datasets. This sort of data is inherently noisy, mainly due to the data acquisition process (drifts, imprecision and shading, etc.) and sensor noise. Furthermore, due to the nature of the capture devices and in particular to the number of sensors that are operated, the nature of the reconstructed skeletons is also subject to some variation for two main reasons: in the one hand the morphology of the actor (segment lengths) is the main source of variability for a given capture device, and in the second hand, the number of joints (the number of degree of freedom) vary according to the capture device, leading to some additional noise due to the reconstruction of the skeletal data from sensor outputs. Thus, any movement can be defined as a multivariate state vector describing a trajectory over time, i.e. a time series

{Y t ∈ R k } T 1 = [Y 1 , ..., Y T ],
where the k spatial dimension (k = 3 • N , with N the number of joints) typically varies between 20 and 100 according to the capture devices and the considered task. As this state vector is obviously not composed of independent scalar dimensions, the spatio-temporal encoded redundancies open new prospects for dimension reduction as well as noise reduction approaches, which is particularly relevant when considering motion recognition as one can targets gains both in terms of computation time and error rate.

III. RELATED WORKS

Gesture analysis and recognition is very broad and recently very active area due to the democratization of low cost motion capture systems by camera. It covers aspects of signal and image processing, dynamic modeling, statistical and machine learning approaches. We hereinafter give a brief non-exhaustive overview of the main methods proposed for gesture analysis and recognition. These methods primarily focus on the extraction of features that significantly represent the kinematics and dynamics of the skeleton data characterized as a whole, or portions of it. Among them, some approaches, based on linear dynamic models [START_REF] Veeraraghavan | Role of shape and kinematics in human movement analysis[END_REF], have exploited autoregressive (AR) models and autoregressive moving-average (ARMA) models to characterize the kinematics of movements, while other approaches, based on nonlinear dynamic models [START_REF] Bissacco | Classification and recognition of dynamical models: The role of phase, independent components, kernels and optimal transport[END_REF], have developed movement analysis and recognition scheme based on dynamical models controlled by Gaussian processes. [START_REF] Mitra | Gesture recognition: A survey[END_REF] propose a synthesis of the major gesture recognition approaches relying on Hidden Markov Models (HMM). [START_REF] Wang | Hidden conditional random fields for gesture recognition[END_REF] have exploited conditional random fields to model joint dependencies and thus increase the discrimination of HMM-like models. Recurrent neural network models have also been used [START_REF] Martens | Learning recurrent neural networks with hessian-free optimization[END_REF]; among them, conditional restricted Boltzman's machines [START_REF] Larochelle | Learning algorithms for the classification restricted boltzmann machine[END_REF] have been studied in the context of motion captured data modeling and classification. Some methods address more specifically the problem of dimensionality reduction with an objective for reducing the variability while seeking an efficiency gain. In particular, the Principal Component Analysis (PCA) has been widely exploited for gesture analysis and recognition with the objective of reducing the dimensionality of motion data [START_REF] Masoud | A method for human action recognition[END_REF]. Other methods such as linear projections preserving locally neighborhoods (Locality Preserving Projection) [START_REF] He | Locality preserving projections[END_REF], or their non-linear counter-parts such as ISOMAP [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF], have been implemented to embed postures in low dimensional spaces in which a more efficient time warp (DTW, see section IV-B) algorithm, associated with the Hausdorff distance, can be used to classify movements. Other ad hoc methods for dimension reduction are also proving their efficiency: we can mention for example the recent work of [START_REF] Yu | Human action recognition with extremities as semantic posture representation[END_REF] that proposes to only consider the trajectories of the 5 end-extremities of the skeleton (2 feet, 2 hands and the head). Models based on Gaussian processes with latent variables are also largely used, for instance a hierarchical version has been recently exploited for gesture recognition [START_REF] Han | Discriminative human action recognition in the learned hierarchical manifold space[END_REF].

The identification of significant variables maximizing the discrimination of motion classes has also been widely explored. [START_REF] Fothergill | Instructing people for training gestural interactive systems[END_REF] and [START_REF] Zhao | Real-time hand gesture detection and recognition by random forest[END_REF] have in particular applied random forests to recognize actions, using a Kinect sensor, while [START_REF] Ofli | Sequence of the most informative joints (smij): A new representation for human skeletal action recognition[END_REF] recently proposed to automatically select the most informative skeletal joints to explain the current action. In the same line, [START_REF] Hussein | Human action recognition using a temporal hierarchy of covariance descriptors on 3d joint locations[END_REF] consider covariance matrices evaluated on some skeletal joints as discriminative descriptors to characterize a movement sequence. The use of sliding windows can be view as a dimension reduction along the time axis. In [START_REF] Li | Action recognition based on a bag of 3d points[END_REF] a simple bag of 3D points is used to represent and recognize gestural action. Similarly, in [START_REF] Wang | Mining actionlet ensemble for action recognition with depth cameras[END_REF], actionlets are defined from Fourrier's coefficients to characterize the most discriminative joints. Histograms of oriented 4D normals have been also proposed in [START_REF] Oreifej | Hon4d: Histogram of oriented 4d normals for activity recognition from depth sequences[END_REF] for the recognition of gestural actions and movements from sequences of depth images.

Finally, it can be mentioned, among many existing applications that address the use of elastic distances into a recognition process, the recent work described in [START_REF] Sempena | Human action recognition using dynamic time warping[END_REF], as well as the hardware acceleration proposed in [START_REF] Hussain | User independent hand gesture recognition by accelerated dtw[END_REF]. However, to our knowledge, no work exploiting this type of distance has directly studied the question of data reduction along the time axis.

IV. DOWNSAMPLING MOVEMENT SEQUENCES COUPLED TO ELASTIC KERNEL MACHINES

When considering the use of elastic distances or kernels to benefit from their ability to deal with some form of temporal variability, we are rapidly confronted with their computational cost, in general quadratic with the length of the time series that are processed and linear with the spatial dimension (number of degrees-of-freedom). This high computational complexity is somehow limiting their use, especially when large amounts of data has to be processed, or when so-called real-time constraint is required. It is therefore a priori particularly relevant to consider a dual dimensionality reduction, firstly on the time axis, and secondly on the spatial axis. Hence, in the context of mocap data processing, it seems useful to determine if a spatio-temporal redundancy of motion paths can be exploited, especially in the perspective of using elastic distances.

A. Dimension reduction along the time axis

Considering the quite rich literature on gesture recognition, it is significant to note that while some studies have shown success with dimensionality reduction on the spatial axis, very few have directly addressed a reduction in dimensionality along the time axis per se to reduce the complexity of elastic matching. [START_REF] Keogh | Scaling up dynamic time warping for datamining applications[END_REF] explicitly mentioned a temporal sub-sampling associated with a dynamic time warping in the context of time series mining, followed later by [START_REF] Marteau | Time warp edit distance with stiffness adjustment for time series matching[END_REF]. In order to explicitly reduce dimensionality along the time axis, our straightforward approach here consists in sub-sampling the motion data so that each motion trajectory takes the form of a fixed-size sequence of L skeletal postures, evenly distributed along the time axis. It becomes then easy to perform a classification or recognition task by using elastic kernel machines on such fixed-size sequences to assess performance rates depending on the degree of sub-sampling that is considered. Indeed, this approach is quite raw, as long sequences can be characterized with the same number of skeletal poses than short sequences. For very short sequences, whose length is shorter than L, if any, we over-sample the sequence in order to meet the fixedsize requirement. But we consider this case as very marginal here since we seek a sub-sampling rate much lower than the average length of the motion sequence.

B. Elastic kernels and their regularization

Dynamic Time Warping (DTW), [START_REF] Velichko | Automatic recognition of 200 words[END_REF], [START_REF] Sakoe | A dynamic programming approach to continuous speech recognition[END_REF], by far the most used elastic measure, is defined as

d dtw (X p , Y q ) = d 2 E (x(p), y(q)) (1)
+ Min

d dtw (X p-1 , Y q ) sup d dtw (X p-1 , Y q-1 ) sub d dtw (X p , Y q-1 )
ins

where d E (x(p), y(q) is the Euclidean distance (possibly the square of the Euclidean distance) defined on R k between the two postures in sequences X and Y taken at times p and q respectively. Besides the fact that this measure does not respect the triangle inequality, it does not directly define a positive definite kernel. When performed by a support vector machine (SVM) model, the optimization problem inherent to this type of learning algorithm is no longer quadratic. Moreover, the convergence towards the optimorum is no longer guaranteed, which, depending on the complexity of the task may be considered as detrimental. Besides the fact that the DTW measure does not respect the triangle inequality, it is furthermore not possible to directly define a positive definite kernel from it. Hence, the optimization problem, inherent to the learning of a kernel machine, is no longer quadratic which could, at least on some tasks, be a source of limitation.

Regularized DTW: recent works [START_REF] Cuturi | A Kernel for Time Series Based on Global Alignments[END_REF], [START_REF] Marteau | On constructing positive elastic kernels with application to time series classification[END_REF] allowed to propose new guidelines to regularize kernels constructed from elastic measures such as DTW. A simple instance of such regularized kernel, derived from [START_REF] Marteau | On constructing positive elastic kernels with application to time series classification[END_REF] for time series of equal length, takes the following form, which relies on two recursive terms :

K rdtw (X p , Y q ) = K xy rdtw (X p , Y q ) + K xx rdtw (X p , Y q ) K xy rdtw (X p , Y q ) = 1 3 e -νd 2 E (x(p),y(q))    h(p -1, q)K xy rdtw (X p-1 , Y q ) h(p -1, q -1)K xy rdtw (X p-1 , Y q-1 ) h(p, q -1)K xy rdtw (X p , Y q-1 ) K xx rdtw (X p , Y q ) = 1 3    h(p -1, q)K xx rdtw (X p-1 , Y q )e -νd 2 E (x(p),y(p)) ∆ p,q h(p, q)K xx rdtw (X p-1 , Y q-1 )e -νd 2 E (x(p),y(q)) h(p, q -1)K xx rdtw (X p , Y q-1 )e -νd 2 E (x(q),y(q)) (2)
where ∆ p,q is the Kronecker's symbol, ν ∈ R + is a stiffness parameter which weights the local contributions, i.e. the distances between locally aligned positions, and d E (., .) is a distance defined on R k .

The initialization is simply

K xy rdtw (X 0 , Y 0 ) = K xx rdtw (X 0 , Y 0 ) = 1.
The main idea behind this line of regularization is to replace the operators min and max (which prevent the symmetrization of the kernel) by a summation operator ( ). This leads to consider, not only the best possible alignment, but also all the best (or nearly the best) paths by summing up their overall cost. The parameter ν is used to control what we call nearly-the-best alignment, thus penalizing more or less alignments too far from the optimal ones. This parameter can be easily optimized through a cross-validation.

Elastic kernels: we consider in this paper only the exponential kernel (Gaussian or RBF-type) constructed from the two previous elastic measures d dtw and K rdtw , and the non-elastic kernel obtained from the Euclidean distance1 , i.e. K dtw (., .) = e -d dtw (.,.)/σ , and K E (., .) = e -d 2 E (.,.)/σ . For the regularized DTW kernel, a data dependent normalization heuristic is required and the final kernel takes the form: K rdtw (., .) = e βK α rdtw (.,.)/σ , with

• α = 1/log(max(K rdtw (., .))/min(K rdtw (., .))) and

• β = exp(-α • log(min(K rdtw (., .)))),
where min and max are taken over all the training data pairs.

V. EXPERIMENTATION

To estimate the robustness of the proposed approach, we evaluate it on two motion capture databases of opposite quality, the first one developed at the Max Planck Institute, the other at Microsoft research laboratories. HDM05 data set [START_REF] Müller | Documentation mocap database hdm05[END_REF] consists of data captured at 120hz by a Vicon MX system composed of a set of reflective optical markers followed by six high-definition cameras and configured to record data at 120hz. The movement sequences are segmented and transformed into sequences of skeletal poses consisting of N = 31 joints, each associated to a 3D position (x, y, z). In practice the position of the root of the skeleton (located near its center of mass) and its orientation serving as referential coordinates, only the relative positions of the remaining 30 joints are used, which leads to represent each position by a vector Y T ∈ R k , with k = 90. We consider two recognition/classification tasks: HDM05-1 and HDM05-2 that are respectively those proposed in [START_REF] Ofli | Sequence of the most informative joints (smij): A new representation for human skeletal action recognition[END_REF] (also exploited in the work of [START_REF] Hussein | Human action recognition using a temporal hierarchy of covariance descriptors on 3d joint locations[END_REF]) and [START_REF] Ofli | Sequence of the most informative joints (smij): A new representation for human skeletal action recognition[END_REF]. For both tasks, three subjects are involved during learning and two separate subjects are involved during testing. For task HDM05-1, 11 gestural actions are processed: {deposit floor, elbow to knee, grab high, hop both legs, jog, kick forward, lie down floor, rotate both arms backward, sneak, squat, and throw basketball}. This constitutes 249 motion sequences. For task HDM05-2 , the subjects are the same, but five additional gestural actions are considered in addition to the previous 11: {jump, jumping jacks, throw, sit down, and stand up}. For this task, the data set includes 393 movement sequences in total. For both tests, the lengths of the gestural sequences are between 56 and 901 postures (corresponding to a movement duration between 0.5-7.5 sec.) .

MSR-Action3D data set:

This database [START_REF] Li | Action recognition based on a bag of 3d points[END_REF] has recently been developed to provide a Kinect data benchmark. It consists of 3D depth image sequences (depth map) captured by the Microsoft Kinect sensor. It contains 20 typical interaction gestures with a game console that are labeled as follows high arm wave, horizontal arm wave, hammer, hand catch, forward punch, high throw, draw x, draw tick, draw circle, hand clap, two hand wave, side-boxing, bend, forward kick, side kick, jogging, tennis swing, tennis serve, golf swing, pickup & throw . Each action was carried out by 10 subjects facing the camera, 2 or 3 times. This data set includes 567 motion sequences whose lengths vary from 14 to 76 skeletal poses. The 3D images of size 640 × 480 were captured at a frequency of 15hz. From each 3D image a skeletal posture has been extracted with N = 20 joints, each one being characterized by three coordinates. As for the previous data set, we characterize postures relatively to the referential coordinates located at the root of the skeleton, which leads to represent each posture by a vector Y t ∈ R k , with k = 3×19 = 57. The task is to provide a cross-validation on the subjects, i.e. 5 subjects participating in learning and 5 subjects participating in testing, considering all possible configurations which represent 252 learning/testing pairs in total.

A. Results and analysis

For the two considered tasks, we present the results obtained using a SVM classifier built from the LIBSVM library [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF], the elastic kernels K dtw and K rdtw , and as a baseline the Euclidean distance kernel, K E . Fig. 2. Classification accuracies for HDM05-1 task (top), and HDM05-2 task (bottom), when the number of skeletal poses varies: K E (red, circle, dash), K dtw (black, square, plain), K rdtw (blue, star, dotted).

Figure 2 presents the classification accuracies when the number of skeletal postures selected after downsampling varies between 5 and 30. Results for HDM05-1 task is presented in the top sub-figure, while results for HDM05-2 task is given in the bottom sub-figure. Figure 3 shows also the classification accuracies obtained on the MSRAction3D data set when the number of skeletal postures varies between 10 and 30: the accuracies obtained for a 10-fold cross validation on the training data is given in the top sub-figure. accuracies for the testing data is given in the bottom sub-figure.

On both figures, we observe that the sub-sampling does not catastrophically degrade the accuracies. High levels of down-sampling (e.g. 10 to 15 postures retained by movement, which represents an average compression ratio of 97 % for HDM05 and 70 % on MSRAction3D) lead to very satisfactory results (96-98 % for the two HDM05 tasks and 95 to 97 % for the MSRAction3D task on the learning data). The SVM classifier constructed on the basis of the regularized kernel K rdtw produces the best recognition rate. We note that the MSRAction3D task is more difficult: much lower performance are obtained for the SVM built on the basis of the Euclidean distance; in addition, if very good classification rate (96 %) is obtained on the training data, due to the noisy nature of Kinect data and the inter subject variability, the recognition rate on the test data drop down to 82 % . For comparison, table II gives results obtained by different methods of the state-of-the-art and compare them with the performance of our SVM constructed from the Regularized DTW associated with a down-sampling of 15 postures. To that end, we have reimplemented the Cov3DJ approach [START_REF] Hussein | Human action recognition using a temporal hierarchy of covariance descriptors on 3d joint locations[END_REF] to get, for the MSRAction3D data set, the average result given by a 5-5 cross-validation on the subjects (252 tests). This comparative analysis shows that the SVM constructed from regularized DTW kernel provides results slightly above the current stateof-the-art for the considered data sets and tasks.

HDM05-1

Accuracy (%) SMIJ [START_REF] Ofli | Sequence of the most informative joints (smij): A new representation for human skeletal action recognition[END_REF] 84.40 Cov3DJ, L = 3 [START_REF] Hussein | Human action recognition using a temporal hierarchy of covariance descriptors on 3d joint locations[END_REF] 95 Finally, in Figure 4, we give the average CPU elapsed time for the processing of a single gestural MSRAction3D action when varying the number of retained skeletal poses. The test has been performed on an Intel Core i7-4800MQ CPU, 2.70GHz. Although the computational cost for the elastic kernel is quadratic, the latency for the classification of a single gestural action is less than 25 milliseconds when 15 poses are considered, which effectively meets easily real-time requirements.

VI. CONCLUSION AND PERSPECTIVES

In the context of isolated gesture recognition, where few studies explicitly consider dimension reduction along the time axis, we have presented a simple approach based on subsampling motion sequences coupled to the exploitation of elastic kernel machines. On the data sets and tasks that we have addressed, we have shown that, even when quite important down-sampling is considered, the recognition accuracy only slightly degrades. The temporal redundancy is therefore high and apparently not critical for the discrimination of the selected movements and tasks. In return, the down-sampling benefits in terms of computational complexity is quadratic with the reduction of the number of skeletal postures kept along the time axis. Furthermore, the elasticity of the kernel provides a performance gain (compared to kernel based on the Euclidean distance) which is important when the data are characterized by high variability. Our results show that a SVM based on a regularized DTW kernel is very competitive comparatively to the state-of-the-art methods applied on the two tested data sets, even when the dimension reduction on the time axis is important. This study opens perspectives to the use of more sophisticated elastic kernels [START_REF] Marteau | Time warp edit distance with stiffness adjustment for time series matching[END_REF] associated to adaptive sampling techniques [34] [35] capable of extracting the most significant and discriminant skeletal poses in movement sequences.

Fig. 1 .

 1 Fig. 1. Examples of skeletons reconstructed from motion data captured from the Kinect (left) and from the Vicon-MX device used by the Max Planck Institute (right).

Fig. 4 .

 4 Fig. 4. Elapsed time as a function of the number of skeletal poses (10 to 30 poses): i) Euclidean Kernel, Red/round/dotted line, ii) Elastic kernel (RDTW), Black/square/plain line.

Table I

 I gives for the MSRAction3D data set and for the

		KE A	KE T	K dtw A	K dtw T	K rdtw A	K rdtw T
	Mean	87,71	69,73	96,04	81,41	96,65	82,50
	Stand. dev.	2,34	5,73	1,36	5,04	1,13	3,22

TABLE I .

 I MEANS AND STANDARD DEVIATIONS OF CLASSIFICATION ACCURACIES ON THE MSRACTION3D DATA SET OBTAINED ACCORDING TO A CROSS-VALIDATION ON THE SUBJECTS (252 TESTS) A: ON THE TRAINING DATA, T: ON TEST DATA FOR A NUMBER OF SKELETAL POSTURES EQUAL TO 15.

		98
		96
	Recognition accuracies (%)	90 92 94
		88
		86
		10	12	14	16	18	20	22	24	26	28	30
							# Skeletal poses
		86
		84
		82
	Recognition accuracies (%)	74 76 78 80
		72
		70
		68
		10	12	14	16	18	20	22	24	26	28	30
							# Skeletal poses
	Fig. 3. Classification accuracies for the MSRAction3D data set, on the
	learning data (top) and testing data (bottom) when the number of skeletal
	poses varies: K means and
	standard deviations, obtained on the training data (A) and
	testing data (T), of recognition rates (classification accuracies)
	when performing the cross-validation over the 10 subjects (252
	configurations). For this test, movements are represented as
	sequences of 15 skeletal postures. The drop of accuracies
	between Learning and testing is due, on this dataset, to the
	large inter subjects variability of movement realizations.

E (red, circle, dash), K dtw (black, square, plain), K rdtw (blue, star, dotted). SVM based on K E , K dtw and K rdtw kernels,

TABLE II .

 II COMPARATIVE STUDY.

	10	12	14	16	18	20	22	24	26	28	30

The Euclidean distance is usable only because a fixed number of skeletal positions is considered to characterize each movement, and this, irrespectively of their initial length


V Jaouen

P Gonzalez

D Guilloteau

I Buvat

C Tauber

4DGVF-BASED FILTERING OF VECTOR-VALUED IMAGES

Keywords: PDE, Image restoration, Vector-valued images, Gradient vector flow, Positron emission tomography

In this paper, we propose a new method for vector-valued image restoration in order to reduce noise while simultaneously sharpening vector edges. Our approach is a coupled anisotropic diffusion and shock filtering scheme that exploits a new robust 4DGVF vector field tailored for vector-valued images. The proposed scheme sharpens edges in directions diffused from the entire spatio-spectral information available with a more precise and a more stable sharpening effect along the iterative processing. We validate our method on color images as well as on realistic simulations of dynamic PET images.

INTRODUCTION

Vector-valued images occur in various contexts such as color images, hyperspectral images (e.g. radar, LiDAR), or medical images acquired at different time intervals (e.g. dynamic positron emission tomography (PET) imaging or functional magnetic resonance imaging). Such images require dedicated methods that take profit of additional information provided by the extra dimension available [START_REF] Sapiro | Anisotropic diffusion of multivalued images with applications to color filtering[END_REF]. In particular, edgebased techniques extended to the vector-valued case typically require a redefinition of the notion of gradient. A popular generalization of this notion to vector-valued images or vector gradient, was proposed by Di Zenzo, following considerations from Riemannian geometry [START_REF] Zenzo | A note on the gradient of a multi-image[END_REF]. According to Di Zenzo, the image I is considered as a vector field whose components are the different channels. The gradient direction is then associated with the direction that maximizes a quadratic form of the total differential dI of I. This idea was later developed using the tensor formalism [START_REF] Lee | Detecting boundaries in a vector field[END_REF], converting this geometric problem into an algebraic one in which the gradient direction can be associated to the dominant eigenvector of a local structure tensor, or second-moment matrix of the image. The gra-

The research leading to these results has received funding from the European Union's Seventh Framework Programme (FP7/2007-2013) under grant agreement HEALTH-F2-2011-278850 (INMiND). dient magnitude is then generally obtained from a combination of the different eigenvalues of the structure tensor. Over the years, Di Zenzo's gradient and the local structure tensor formalism have been used for different applications such as feature or saliency detection [START_REF] Van De Weijer | Robust photometric invariant features from the color tensor[END_REF][START_REF] Piella | Image fusion for enhanced visualization : A variational approach[END_REF], image segmentation [START_REF]Vector (self) snakes : A geometric framework for color, texture, and multiscale image segmentation[END_REF][START_REF] Xie | RAGS : Region-aided geometric snake[END_REF], edge detection (known as vector edge) [START_REF] Cumani | Edge detection in multispectral images[END_REF], and image filtering [START_REF] Sapiro | Anisotropic diffusion of multivalued images with applications to color filtering[END_REF][START_REF] Chan | Active contours without edges for vector-valued images[END_REF][START_REF] Tschumperlé | Diffusion PDEs on vector-valued images[END_REF][START_REF] Tschumperlé | Vector-valued image regularization with PDEs : A common framework for different applications[END_REF].

For single channel images, many non-linear filtering approaches, such as the anisotropic diffusion filter of Perona and Malik (PM) and its refinements [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF][START_REF] Catté | Image selective smoothing and edge detection by nonlinear diffusion[END_REF][START_REF] Kornprobst | Nonlinear operators in image restoration[END_REF] have been proposed in order to reduce noise while preserving edge strength. Alvarez and Mazzora proposed a filtering equation that couples PM anisotropic diffusion and a shock filter term [START_REF] Alvarez | Signal and image restoration using shock filters and anisotropic diffusion[END_REF][START_REF] Osher | Feature-oriented image enhancement using shock filters[END_REF] in order to further enhance the edge signal and restore sharpness. This term can be seen as an inverse diffusion that restores abrupt discontinuities around edges. However, a direct extension of this latter approach to the vector-valued case, by applying it independently on each channel, is not satisfying. Indeed, edges are often difficult to identify accurately in individual channels of vector-valued images, and independent morphological operations acting in the different channels are likely to create false spectral characteristics, such as false colors [START_REF] Serra | The "false colour" problem[END_REF]. Based on the geometrical considerations of Di Zenzo, Tschumperlé and Deriche extended the Alvarez-Mazorra formulation to the vector-valued case by combining a nonlinear diffusion term, a shock filter term and a data-fidelity term for color image restoration [START_REF] Tschumperlé | Diffusion PDEs on vector-valued images[END_REF]. By acting on all channels in the same fashion, this method reduces incoherent restoration of the vector components. However, shock filters rely on the localization of the Laplacian zerocrossings and therefore remain very sensitive to noise, even when coupled to vector diffusion schemes.

In this paper, we propose a novel approach for vectorvalued image sharpening and denoising based on a coupled vector diffusion and shock filtering scheme. The proposed approach exploits sharpening directions drawn from a new 4D gradient vector flow (4DGVF) field that propagates and regularizes vector gradient information throughout the image. By partially transferring the vector edges localization and orientation tasks the 4DGVF field, we define more accurate and robust directions along which blur is compensated in a stable fashion. As all image channels are sharpened in coherent directions, the proposed method also reduces the appearance of spurious spectral characteristics. To a certain extent, this approach can be considered as a generalization of the GVFbased anisotropic diffusion model proposed by Yu and Chua for 2D grayscale images to the vector-valued case [START_REF] Yu | GVF-based anisotropic diffusion models[END_REF]. We validate our method on synthetic color images and realistic simulations of 4D (3D+t) dynamic PET images, and compare our results to the approaches of Yu and Chua [START_REF] Yu | GVF-based anisotropic diffusion models[END_REF] and of Tschumperlé and Deriche [START_REF] Tschumperlé | Diffusion PDEs on vector-valued images[END_REF].

PROPOSED METHOD

In this section, we present our vector-valued restoration approach based on the 4DGVF field, a gradient vector flow (GVF) field tailored for vector-valued images [START_REF] Xu | Generalized gradient vector flow external forces for active contours[END_REF]. We start by defining a vector geometry in order to characterize vector edges at a local scale, both in amplitude and in direction. These geometrical considerations are the foundation upon which the 4DGVF field and the proposed regularization scheme are built.

Definition of a vector geometry

In the continuous domain, we denote by I a nD vectorvalued image constituted of M channels :

I(x, c) : (Ω s ⊗ Ω c ) ∈ R n ⊗ N → R,
where Ω s is the nD spatial domain of the image and Ω c the channel dimension. x = (x 1 , ..., x n ) ∈ Ω s is the voxel position in each channel. We denote by I k the k th channel of the image.

According to Di Zenzo, I is a nD→ M D vector field. The gradient field of the multichannel image is locally oriented in the direction which maximizes the quadratic form of the total differential dI of I [START_REF] Zenzo | A note on the gradient of a multi-image[END_REF] :

dI = n i=1 ∂I ∂x i dx i . (1)
We define a weighted quadratic form for dI, or first fundamental form :

dI 2 ω = dx T G ω dx, (2)
with G ω a regularized, weighted structure tensor of the image [START_REF] Weickert | Anisotropic diffusion in image processing[END_REF] G

ω = K σ * M k=1 ω k (∇I k ⊗ ∇I k T ), (3)
where K σ is a Gaussian kernel of scale σ, * is the convolution operation, ⊗ is the tensor product, and ω k is a weighting factor for channel I k . The weighting factors allow to characterize the reliability of the different channels. Without a priori estimation, a natural generalization of Di Zenzo's approach is to weight all channels equally by choosing ω k = 1/M, ∀k.

Let λ 1 ≥ ... ≥ λ n be the ordered set of eigenvalues of the structure tensor G ω and θ 1 , ..., θ n the associated orthonormal eigenvector set. The eigenvalues λ i give the scalar rates of change of the first fundamental form in a local basis of extremal variations. The dominant eigenvector θ 1 associated with λ 1 is collinear to the vector gradient, while the remaining eigenvectors span the hyperplane tangent to the local isophote. In order to measure the vector gradient amplitude, several combinations of the eigenvalues can be found in the literature [START_REF] Tschumperlé | Diffusion PDEs on vector-valued images[END_REF]. In our study, we choose a coherence norm proposed by Weickert in [START_REF] Weickert | Coherence-enhancing diffusion of colour images[END_REF] :

dI 2 ω = n-1 i=1 m j=i+1 (λ i -λ j ) 2 . ( 4
)
This norm measures the amount of local anisotropy, a generalization to the nD case of the norm presented by Sapiro in [START_REF] Sapiro | Anisotropic diffusion of multivalued images with applications to color filtering[END_REF]. This measure exhibits oriented gradient patterns in the image. We thus define a measure N ω of the vector edge amplitude.

N ω := dI ω . (5)

4DGVF field

We have recently proposed the 4DGVF vector field for the robust segmentation of vector-valued images with active surfaces [START_REF] Jaouen | 4-D gradient vector flow : segmentation par surface active pour images multi-composantes[END_REF]. This field is obtained through a gradient vector flow-like scheme that propagates vector gradients in the image in a nonlinear fashion [START_REF] Xu | Generalized gradient vector flow external forces for active contours[END_REF]. It exploits both amplitude and directional information contained in the analysis of the structure tensor, on the contrary of previous approaches that rely on the norm of the tensor only [START_REF] Yang | Unsupervised segmentation based on robust estimation and color active contour models[END_REF][START_REF] Jaouen | Vector-based active surfaces for segmentation of dynamic PET images[END_REF].

We define a vectorial edge map V, a vector field collinear to the local dominant eigenvectors of G ω , but oriented toward the nearest vector edges :

V = θ + sign < θ + , ∇N ω > . (6)
The 4DGVF field is obtained by nonlinear diffusion of the components of V throughout the image. It is defined as the steady-state solution F ext of the following vector partial differential equation :

∂ F ext ∂t = g(N ω )∆ F ext -h(N ω )( F ext -V), (7)
where g(s) = e -|∇s|/κ and h = 1 -g are two functions that control the tradeoff between the first and second terms through parameter κ [START_REF] Xu | Generalized gradient vector flow external forces for active contours[END_REF], and ∆ is the vector Laplace operator. In the vicinity of vector edges, as measured by N ω , the directions of F ext are constrained by V, while isotropic diffusion of F ext prevails in homogeneous regions. We obtain a regularized vector field oriented locally toward the nearest vector edge.

4DRSF vector regularization

We propose a new regularization approach for nD vectorvalued images that combines anisotropic diffusion and 4DGVFregularized shock filtering. The proposed 4D regularized shock filter (4DRSF) approach is driven by the following set of coupled partial differential equations :

∂I k ∂t = c D+ (N ω ) ∂ 2 I k ∂ θ 2 1 + c D-(N ω ) n i=2 ∂ 2 I k ∂ θ 2 i + α( 0 I k -I k ) -c S (N ω )sign(< F, ∇I k |∇I k | >) ∂I k ∂ θ 1 , (8)
where operators ∂/∂ θ i and ∂ 2 /∂ θ i 2 respectively correspond to the directional derivative of first and second order in the θ i direction and where < , > is the dot product.

The first two terms in equation ( 8) respectively correspond to a diffusion along the vector gradient direction and along the hyperplane tangent to the local isophote, weighted by diffusion coefficients c D+ and c D-, decreasing functions of N ω . The third term is a classical data-attachment term controlled by parameter α, where 0 I k denotes the original k th channel. The last term in eq. ( 8) is the 4DGVF-based shock filter term. It sharpens vector edges in robust directions established prior to diffusion by the 4DGVF field F ext . c S is a decreasing function of N ω so that sharpening occurs more near edges than in flat regions. While previous shock filter based regularization schemes ( [START_REF] Alvarez | Signal and image restoration using shock filters and anisotropic diffusion[END_REF][START_REF] Tschumperlé | Diffusion PDEs on vector-valued images[END_REF]) sharpen edges in the direction of the curvature zero-crossings, which may change along the iterations, such a scheme guarantees a stable edge-enhancement behavior along robust directions determined by the 4DGVF field.

Implementation

We implemented the 4DRSF regularization method using MATLAB. We used the following diffusion and shock coefficients :

c -(s) = 1 √ 1+ s r , c + = c 2 -and c S = 1-c -,
where r is a scale parameter. In order to handle data with similar orders of magnitude, all intensities were scaled to the range [0 1].

RESULTS

We tested our method on two different data sets for which ground truths were available : synthetic color images and realistic simulations of dynamic PET images.

Comparative evaluation

We compared the 4DRSF approach to three other methods : 1) an isotropic Gaussian filter applied in each channel I k independently ; 2) the GVF-based anisotropic diffusion model of Yu and Chua for single-channel images [START_REF] Yu | GVF-based anisotropic diffusion models[END_REF]. This method couples anisotropic diffusion with a shock filter term that exploits directions established prior to diffusion by a classical GVF vector field. As for the Gaussian filter, this method was applied in each channel independently ; 3) the regularization approach of Tschumperlé and Deriche, which takes profit of Di Zenzo's vector gradient in order to reduce noise and to sharpen edges in a coherent way along the channels of vectorvalued images [START_REF] Tschumperlé | Diffusion PDEs on vector-valued images[END_REF].

The parameters of all methods were established so as to minimize the root mean square error (RMSE) between the result and the ground truth. We evaluated two quantitative criteria : the RMSE and the signal-to-noise ratio (SNR) [START_REF] Pratt | Digital Image Processing : PIKS Inside[END_REF].

Synthetic color images

We generated a dataset of 50 synthetic 2D color images having various levels of Gaussian blur and of additive Gaussian noise. The filtering process was performed in the RGB color space. The ground truth is displayed in Fig. 1a anda representative image of the dataset is shown in Fig. 1b. Results of the different approaches are shown in figures 1c-f. For this image, the isotropic Gaussian filter (Fig. 1c) reduced noise at the expense of lowering the spatial resolution. The anisotropic diffusion filter of Yu and Chua (Fig. 1d) reduced noise rather efficiently across the image while providing a sharp rendering of edges. However, a spurious edge sharpening behavior occurred in flat regions. Moreover, false colors were produced around edges. This can be explained by the fact that edges are not sharpened with the same amplitude or along the same directions in the R, G and B channels. The approach of Tschumperlé and Deriche (Fig. 1e) led to a well-rendered image, mostly avoiding the above mentioned issues. However, vector edges were less sharpened that with the method of Yu and Chua, as shocks are developed in varying directions along the iterative treatment, causing the spatial smearing of edges. A visual inspection suggests that the proposed 4DRSF approach was better able to restore the original piecewise constant aspect of the image. These observations were confirmed by the quantitative results displayed in table 1. For this dataset, the Gaussian filter (Gaus.), the approach of Yu and Chua (Yu.), the approach of Tschumperlé and Deriche (Tsch.) and the proposed 4DRSF approach increased the SNR of the unfiltered image (U.I.) by 42%, 63%, 61% and 65% respectively.

Dynamic PET image simulations

The validation of filtering results using clinical images is difficult due to lack of ground truth. We assessed the 4DRSF regularization approach with quantitative results on realistic Monte Carlo simulations of dynamic PET images. Three images with different levels of noise were generated using GATE, a highly realistic PET image simulator based on the CERN's GEANT4 particle interaction platform [START_REF] Jan | GATE v6 : a major enhancement of the gate simulation platform enabling modelling of CT and radiotherapy[END_REF]. Noise While the isotropic Gaussian filter (Fig. 2c) lowered the spatial resolution, the channel-by-channel approach of Yu and Chua (Fig. 2d) led to better denoising. However, it failed at providing sharp boundaries and the image remained blurred. Indeed, this method does not take advantage of the redundancy of the edge signal along the channels in order to lower its sensitivity to noise. The GVF fields calculated in individual channels were too much affected by noise and thus failed at providing relevant sharpening directions. The vector-valued restoration approach of Tschumperlé and Deriche (Fig. 2e) sharpened edges in a coherent way and benefited of the spectral information to identify and to locate edges more precisely. However, the geometry of the different tissues was altered during the process. For example, the thalamus and the frontal lobe regions have been artificially joined due to spurious diffusion behavior. This can be explained by the displacements of the Laplacian zero-crossings along the iterations. Once again, the 4DRSF approach provided both better denoising and better sharpening ability, a property which was confirmed by the quantitative results displayed in the last column of table 1.

To illustrate the method in a preclinical context, an example on a real [18F]-DPA-714 dynamic PET image is shown in Fig. 3. Quinoleic acid was injected in the right striatum of a rat, producing a hypersignal in the image displayed. The filtering result was consistent with the expected morphology of the pathological region.

CONCLUSION

We have proposed a new restoration method for vectorvalued images that allows to increase the SNR while simultaneously sharpening vector edges. The proposed approach exploits robust sharpening directions derived from a new vector field that takes profit of the entire spatio-spectral information available in order to better identify both edges strength and direction. Results obtained on synthetic images confirm the potentiality of the proposed method for the restoration of 2D color and 3D+t images.

Fig. 1 .

 1 Fig. 1. Results for the synthetic image. (a) Ground truth, (b) unprocessed image, (c) isotropic Gaussian smoothing, (d) Yu and Chua channel-by-channel approach, (e) Tschumperlé and Deriche approach, (f) proposed 4DRSF approach.

Fig. 2 .

 2 Fig. 2. Sagittal slice of a representative channel of a dynamic PET image simulation. (a) Ground truth, (b) unprocessed image, (c) isotropic Gaussian smoothing, (d) Yu and Chua channel-by-channel approach, (e) Tschumperlé and Deriche approach, (f) proposed 4DRSF approach.

Fig. 3 .

 3 Fig. 3. Real [18F]DPA-714 PET acquisition of a rat brain. left : original image. middle,right : 4DRSF filtering result

Table 1 .

 1 Average scores for the two datasets

		2D color images	4D PET simulations
		RMSE	SNR (dB)	RMSE SNR (dB)
	U.I.	.15±.04 11.0±1.9 .27±.09 07.9±3.6
	Gaus.	.14±.03 11.5±1.7 .15±.04 10.7±1.0
	Yu.	.11±.01 14.9±0.3 .12±.01 11.2±0.4
	Tsch.	.10±.01 14.9±0.3 .11±.01 11.7±0.8
	4DRSF .09±.01 15.2±0.3 .10±.01 12.9±1.9