Member, IEEE Kenta Kasai

Member, IEEE David Declercq

Senior, IEEE Kohichi Sakaniwa

Fountain Coding via Multiplicatively Repeated Non-Binary LDPC Codes

Keywords: fountain codes, rateless codes, non-binary LDPC codes

We study fountain codes transmitted over the binary-input symmetric-output channel. For channels with small capacity, receivers in fountain coding systems needs to collects many channel outputs to recover information bits. Since a collected channel output yields a check node in the decoding Tanner graph, the channel with small capacity leads to large decoding complexity. In this paper, we introduce a novel fountain coding scheme with non-binary LDPC codes. The decoding complexity of the proposed fountain code does not depend on the channel. Numerical experiments show that the proposed codes exhibit better performance than conventional fountain codes, especially for moderate number of information bits.

I. INTRODUCTION

In this paper, we study fountain codes transmitted over the binary-input symmetric-output channel with capacity C. Fountain codes are a class of erasure-recovering or error-correcting codes which produce limitless sequence of encoded bits from k information bits so that receivers can recover the k information bits from any (1 + ǫ)k/C encoded bits. We denoted overhead by ǫ. Fountain codes encodes k information bits to infinitely long coded bits. The purpose of fountain coding is to realize reliable transmissions under reception of arbitrary (1 + ǫ)k/C received signals among coded bits of length ∞. Even if there is no errors or erasures, this goal can not be achieved by conventional block codes.

The name fountain is after water fountains which endlessly produce water drops to entertain people. Designing fountain codes with small overhead is desirable. LT codes [START_REF] Luby | LT codes[END_REF] and raptor codes [START_REF] Shokrollahi | Raptor codes[END_REF] are fountain codes which achieve vanishing overhead ǫ → 0 in the limit of large information size over the The material in this paper was presented in part at 6th Int. Symp. on Turbo Codes and Related Topics [START_REF] Kasai | Fountain codes with multiplicatively repeated non-binary LDPC codes[END_REF].

channel with C = 1. Such a channel is viewed as the binary erasure channel (BEC) with zero erasure probability in [START_REF] Luby | LT codes[END_REF], [START_REF] Shokrollahi | Raptor codes[END_REF]. By a nice analogy between the BEC and the packet erasure channel (e.g., Internet), fountain codes have been successfully adopted by several industry standards. In such packet erasure channels with many receivers (e.g., 10 6 ), ack-based protocol and re-transmissions can not be afforded due to high server load. Note that if the channel is not heavily impaired, the re-transmission scheme could be a good choice.

In [START_REF] Etesami | Raptor codes on binary memoryless symmetric channels[END_REF], Etesami et al. investigated raptor codes used over the memoryless binary-input outputsymmetric (MBIOS) channels. And they showed that over the binary-input additive white Gaussian noise (BIAWGN) channels with capacity C ≥ 0.49, raptor codes achieve overhead ǫ ≤ 0.08 at bit error rate 10 -7 with information size k = 65536. A raptor code can be viewed as concatenation of an outer high-rate LDPC code and infinitely many single parity-check codes of length d, where d is chosen randomly with probability Ω d for d ≥ 1. In [START_REF] Venkiah | Jointly decoded raptor codes: analysis and design for the BIAWGN channel[END_REF], Venkiah et al. proposed a joint decoding of the concatenated codes and an optimization method for output degree distributions Ω(x) := d≥1 Ω d x d and showed that the optimized codes outperform the conventional ones.

The problems for constructing fountain codes used for general channels with finite inputs are summarized as follows.

• Problem 1: The output degree distribution Ω(x) needs to be optimized for each k. And large check node degree d leads to the large encoding and decoding complexity. • Problem 2: The number of check nodes in the inner codes is given by (1 + ǫ)k/C. This increases as the channel capacity C decreases.

Since check node computation is dominant in decoding, the decoding complexity is high for small C.

• Problem 3: Conventional research has been focused on fountain codes with large information size and vanishing overhead. This leads to large size of memory devices and transmission latency. Good fountain codes with moderate information size is desired. In this paper, we will propose a novel fountain coding scheme which is free of those drawbacks.

In this paper, we consider non-binary LDPC codes defined by sparse parity-check matrices over GF(2 m ) for 2 m > 2. Non-binary LDPC codes were invented by Gallager [START_REF] Gallager | Low Density Parity Check Codes[END_REF] and, Davey and MacKay [START_REF] Davey | Low-density parity check codes over GF(q)[END_REF] found non-binary LDPC codes can outperform binary ones. Non-binary LDPC codes have captured much attention recently due to their decoding performance.

It is known that the irregularity of Tanner graphs helps in improving the decoding performance of binary LDPC codes. On the other hand, it is not the case for the non-binary LDPC codes. Interestingly, the (2, d c )-regular non-binary LDPC codes over GF(2 m ) are empirically known [START_REF] Poulliat | Design of regular (2,dc)-LDPC codes over GF(q) using their binary images[END_REF] as the best performing codes for 2 m ≥ 64, especially for short code length. This suggests that, for designing non-binary LDPC codes, one does not need to optimize the degree distributions of Tanner graphs, since (2, d c )-regular non-binary LDPC codes are best. Furthermore, sparsity of (2, d c )-regular Tanner graph makes efficient decoding possible.

The rest of this paper is organized as follows. Section II defines the proposed codes. Section III describes the decoding algorithm for the proposed codes. In Section IV, we investigate the asymptotic overhead for the BEC and BIAWGN channels by density evolution [START_REF] Richardson | Design of capacity-approaching irregular low-density paritycheck codes[END_REF], [START_REF] Rathi | Density Evolution, Threshold and the Stability Condition for non-binary LDPC Codes[END_REF] and Monte Carlo density evolution, respectively. In Section V, for the BIAWGN channels, we compare the decoding performance of proposed codes and the best known codes for short and moderate code length so far. In Section VI, we give a conclusion.

II. FOUNTAIN CODING WITH MULTIPLICATIVELY REPEATED NON-BINARY LDPC CODES

In this section we explain a new fountain coding scheme. The new coding scheme uses a non-binary LDPC code as a pre-code.

In [START_REF] Kasai | Multiplicatively repeated non-binary LDPC codes[END_REF], the authors introduced low-rate nonbinary codes. The code is a concatenation of (2, 3)-regular non-binary LDPC code and inner multiplicative repetition codes. In general, low-rate LDPC codes have many check nodes and suffer from the high decoding complexity than high rate codes. One of the remarkable feature of the codes proposed in [START_REF] Kasai | Multiplicatively repeated non-binary LDPC codes[END_REF] is that the decoding complexity does not depend on the coding rate. The code exhibits excellent decoding performance for moderate code length and is rate-compatible. We will use the low-rate code [START_REF] Kasai | Multiplicatively repeated non-binary LDPC codes[END_REF] with vanishing rate as a fountain code.

We fix a Galois field GF(2 m ) with a primitive element α and its primitive polynomial π(x). Once the primitive element is fixed, one can represent each symbol in the Galois field as a binary sequence of length m [START_REF] Macwilliams | The Theory of Error-Correcting Codes[END_REF]. For example, with a primitive element α ∈ GF(2 3 ) such that π(α) = α 3 + α + 1 = 0, each symbol is represented as 0 = (0, 0, 0), 1 = (1, 0, 0), α = (0, 1, 0), α 2 = (0, 0, 1), α 3 = (1, 1, 0), α 4 = (0, 1, 1), α 5 = (1, 1, 1) and α 6 = (1, 0, 1). We refer to elements in GF(2 m ) as symbols if m ≥ 2 and bits if m = 1. In this setting, k information bits can be represented as k/m symbol sequence (s 1 , . . . , s k/m ) ∈ GF(2 m ) k/m . We denote the j-th bit of the binary representation of x i ∈ GF(2 m ) by x (j) i . As mentioned in Section I, fountain coding has an application for the packet erasure channel. In the analogous packet scenario, what we consider as a packet is not a symbol s (j) i ∈ GF(2 m ), but a bit s (j) i ∈ GF [START_REF] Luby | LT codes[END_REF]. We denote the data size of packet by D [bits] which is set to very large, e.g., D = 12000 to mitigate the rate loss due to the header overhead in practice. In this paper, we assume D = 1 for simplicity. The number of information packets is given by k/m.

A non-binary LDPC code C over GF(2 m ) is defined by the null space of a sparse M × N paritycheck matrix H = {h c,v } defined over GF(2 m ).

C = {x ∈ GF(2 m ) N | Hx T = 0 ∈ GF(2 m ) M }
The c-th parity-check equation for c = 1, . . . , M is written as

h c,1 x 1 + • • • + h c,N x N = 0 ∈ GF(2 m ), where h c,1 , . . . , h c,N ∈ GF(2 m ) and x 1 , . . . , x N ∈ GF(2 m ).
The binary LDPC codes are represented by Tanner graphs with variable and check nodes [13, pp. 75]. The non-binary LDPC codes, in this paper, are also represented by bipartite graphs with Fig. 1. An example of a pre-code C1. A non-binary (2,3)-regular LDPC code of rate 1/3 over GF(2 m ). Each variable node represents a symbol in GF(2 m ). Each check node represents a parity-check equation over GF (2 m ). The code length is 18 symbols in GF(2 m ) or equivalently 18m bits.

variable nodes and check nodes, which are also referred to as Tanner graphs. For a given sparse parity-check matrix H = {h c,v } over GF(2 m ), the graph is defined as follows. The v-th variable node and c-th check node are connected if h c,v = 0. By v = 1, . . . , N and c = 1, . . . , M , we also denote the v-th variable node and c-th check node, respectively.

A non-binary LDPC code with a parity-check matrix over GF(2 m ) is called (d v , d c )-regular if all the columns and all the rows of the paritycheck matrix have weight d v and d c , respectively, or equivalently all the variable and check nodes have degree d v and d c , respectively. Let C 1 be a (2, 3)-regular LDPC pre-code defined over GF(2 m ) of length N symbols or equivalently mN bits and of rate 1/3. It can be seen that N = 3k/m. The precode C 1 has a 2N/3 × N sparse parity-check matrix H = {h c,v } over GF(2 m ). The matrix H has row weight 3 and column weight 2. Fig. 1 shows the Tanner graph of C 1 of length N =18 symbols. It can be shown that (2, d c )-regular non-binary LDPC codes is linear-encodable by using a non-singular zig-zag sub-graph.

We define a new fountain code C ∞ : GF(2) k → GF(2) ∞ by using the following encoding procedure. 1) First, map k information bits (s

(1) 1 , . . . , s (m) k/m ) to k/m information GF(2 m )-symbols (s 1 , . . . , s k/m ). 2) By the pre-code C 1 , encode the k/m informa- tion symbols to N symbols (x 1 , . . . , x N ) ∈ GF(2 m ) N . 3) Repeat the followings endlessly from i = 1 to ∞. a) Pick uniformly at random v i ∈ [1, N ], w i ∈ [1, m] and r i ∈ GF(2 m ) \ {0}. b) Transmit w i -th bit of r i x v i ∈ GF(2 m ).
The proposed fountain code C ∞ can be viewed as a non-binary raptor code with a non-binary (2,3)-regular LDPC pre-code and an output degree distribution Ω(x) = x [START_REF] Shokrollahi | Raptor codes[END_REF]. Note that Ω(x) = x does not mean simple repetition of bits but multiplicative repetition of symbols in GF(2 m ) for the proposed non-binary setting.

It may be possible to improve the overhead by optimizing general output degree distributions Ω(x) = dmax d=1 x d with some bounded maximum degree d max . In this paper, we focus on the output degree distribution as Ω(x) = x.

III. DECODING SCHEME

We assume that transmission takes place over the MBIOS channel. Specifically, the channel is specified by the transition probability P (•|•) such that P (y|x) = Pr(Y = y|X = x) where X and Y are the random variable of an input bit x and the channel output y, respectively. We assume that the information bits are chosen with uniform probability.

The most important feature of the fountain coding system is that the decoder does not receive all the channel outputs but collects n channel outputs. The decoder recovers the k information bits from the n collected channel outputs. Overhead ǫ is defined [START_REF] Etesami | Raptor codes on binary memoryless symmetric channels[END_REF], [START_REF] Venkiah | Jointly decoded raptor codes: analysis and design for the BIAWGN channel[END_REF] by

ǫ = C/R -1, R = k/n,
where C is the channel capacity. Then, the decoder has n = (1 + ǫ)k/C collected channel outputs.

Note that, in the original setting of fountain codes as in [START_REF] Luby | LT codes[END_REF], [START_REF] Shokrollahi | Raptor codes[END_REF], the capacity is set C = 1, i.e., all the collected bits are uncorrupted. The aim of fountain coding in this paper is to reliably recover the information bits with small overhead. The overhead ǫ = 0 implies that the information bits are transmitted at rate R = C, which is our extreme aim. With infinitely many information bits, raptor codes can achieve ǫ = 0 for the channel with C = 1, i.e., the BEC with zero erasure probability. And raptor codes optimized for the BEC exhibit a quite good performance with large information bits k = 65536. However, for both the BEC and the general MBIOS channels with C < 1, raptor codes exhibit high error floors [START_REF] Etesami | Raptor codes on binary memoryless symmetric channels[END_REF], [START_REF] Venkiah | Rate splitting issue for finite length Raptor codes[END_REF], [START_REF] Venkiah | Analyse et optimisation de codes Raptor pour les canaux multicast sans fil[END_REF], [START_REF] Venkiah | Jointly decoded raptor codes: analysis and design for the BIAWGN channel[END_REF] for moderate information bits with k ≈ 1024.

For the i-th transmitting bit, the sender picked

v i ∈ [1, N ], w i ∈ [1, m] and r i ∈ GF(2 m ) \ {0}
uniformly at random and transmitted w i -th bit of r i x v i ∈ GF(2 m ). This computation r i x v i can be viewed as concatenation with a rate-1 inner random code of length m bits. Let I be the set of Fig. 2. An example of a Tanner graph used for decoding. Some variable nodes are of degree one. The variable nodes of degree one are corresponding to the transmitted symbols whose channel outputs are collected by the decoder. White dots represent bits corresponding to the received channel outputs. Each transmitted bit is encoded in the scheme explained in Section II. It can be seen that the decoder collected 22 channel outputs for this example.

transmitting indices that the receiver collected. It follows #I = n, where #I denotes the cardinality of I. In other words, for i ∈ I, the receiver collects y i that is the corrupted version of the ith transmitted bits. We assume that the decoder knows not only y i but also the indices v i , w i and the multiplicative coefficients r i for i ∈ I. In practice, this is realized by embedding the indices in the header of packets or synchronization between the sender and the receivers [START_REF] Shokrollahi | Raptor codes[END_REF]. If the data size of packet is sufficiently large (e.g., D=12,000 [bits]), the rate loss due to the header is negligible.

The proposed code C ∞ can be decoded by the sum-product (SP) decoding algorithm on the Tanner graphs. The SP decoder for the non-binary LDPC codes exchanges probability vectors in ∈ R 2 m , called messages, between variable nodes and check nodes [START_REF] Rathi | Density Evolution, Threshold and the Stability Condition for non-binary LDPC Codes[END_REF]. An example of the Tanner graph used by the decoder is shown in Fig. 2. The variable nodes of degree one with white dots in Fig. 2 represent collected channel outputs. If the SP decoding algorithm is immediately applied to the proposed codes, all the variable nodes and check nodes, including the variable nodes of those multiplicative repetition symbols, are activated, i.e. exchage the messages. However, the messages reached at the variable nodes of degree one do not change messages that sent back from the nodes. Therefore, after the initialization, the decoder does not need to pass the messages all the way to those variable nodes of degree 1 and their adjacent check nodes of degree 2. Consequently, the decoder uses only the Tanner graph of the pre-code C 1 . It follows that the complexity of the decoding algorithm does not depend on the number n of collected channel outputs and the channel capacity C. In contrast the decoding complexity of the conventional fountain codes largely depends on n and C as explained in Section I.

The SP decoding involves mainly 4 parts, i.e. the initialization, the check to variable computation, the variable to check computation, and the tentative decision parts. Let X be the random variable of a transmitted bit x, and let Y be the random variable of the corresponding channel output y. The a posterior probability Q(x|y) := Pr(X = x|Y = y), for x = 0, 1 and y ∈ A is assumed to be known to the decoder, where A is the receiving alphabet. initialization : The decoders collected n = (1 + ǫ)k/C channel outputs, y i for i ∈ I, where #I = n. Define

I v := {i ∈ I | v i = v}. It follows that I = ∪ N v=1 I v . For each variable node v in C 1 for v = 1, . . . , N , calculate p (0) v (x) for x ∈ GF(2 m ) as follows. p (0) v (x) = ξ i∈Iv p(0) i (r i x) (1)
p(0) i (x) = Q(0|y i ) if the w i -th bit of x is 0 Q(1|y i ) if the w i -th bit of x is 1,
where ξ is the normalized factor such that

(ℓ+1) cv ∈ R 2 m to each adjacent variable node v ∈ ∂c. p(ℓ) vc (x) = p (ℓ) vc (h -1 v,c x) for x ∈ GF(2 m ), p(ℓ+1) cv = ⊗ v ′ ∈∂c\{v} p(ℓ) v ′ c p (ℓ+1) cv (x) = p(ℓ+1) cv (h v,c x) for x ∈ GF(2 m ).
where

p 1 ⊗ p 2 ∈ R 2 m is convolution of p 1 ∈ R 2 m
and p 2 ∈ R 2 m . To be precise,

(p 1 ⊗p 2 )(x) = y,z∈GF(2 m ) x=y+z p 1 (y)p 2 (z) for x ∈ GF(2 m ).
The convolution seems the most complex part of the decoding. The convolutions are efficiently calculated via FFT and IFFT [START_REF] Declercq | Decoding algorithms for nonbinary LDPC codes over GF(q)[END_REF], [START_REF] Rathi | Density Evolution, Threshold and the Stability Condition for non-binary LDPC Codes[END_REF]. Update the iteration round as ℓ := ℓ + 1.

variable to check output : Each variable node v = 1, . . . , N in C 1 has 2 adjacent check nodes since the pre-code C 1 is (2, 3)regular. Let ∂v be the set of adjacent check nodes. The message p

(ℓ) vc ∈ R 2 m sent from v to c ∈ ∂v is given by p (ℓ) vc (x) = p (0) v (x) c ′ ∈∂v\{c} p (ℓ) c ′ v (x) for x ∈ GF(2 m ).

tentative decision

For each v = 1, . . . , N , the tentatively estimated v-th transmitted symbol is given as

x(ℓ) v = argmax x∈GF(2 m ) p (0) v (x) c ′ ∈∂v p (ℓ) c ′ v (x). If x(ℓ) := (x (ℓ) 1 , . . . , x(ℓ) N ) forms a codeword of C 1 , i.e. x(ℓ) satisfies every parity-check equation of C 1 v∈∂c h c,v x(ℓ) v = 0 ∈ GF(2 m )
for all c = 1, . . . , M , the decoder outputs x(ℓ) as the estimated codeword. Otherwise repeat the check to variable, variable to check and tentative steps. If the iteration round ℓ reaches at a pre-determined number, the decoder collects more channel outputs and start over the decoding.

IV. ANALYSIS OF ASYMPTOTIC OVERHEAD

In this section, we investigate the overhead ǫ in the limit of many information bits k → ∞ for the transmissions over the BEC, i.e., C = 1. Rathi developed the density evolution which enables the prediction of the decoding performance of the non-binary LDPC codes in the limit of large code length. The density evolution usually gives, for a given code ensemble, the maximum channel erasure probability, referred to as threshold, at which the average decoding erasure probability goes to zero. We will use the density evolution calculating the smallest overhead ǫ at which the average decoding erasure probability goes to zero in the limit of k → ∞.

The density evolution used in this section was originally developed for the non-binary LDPC code ensembles with parity-check matrices defined over the general linear group GL(GF(2), m). However, Rathi observed that the threshold for the code ensemble defined over GF(2 m ) and GL(GF(2), m) also have approximately the same threshold within the order of 10 -4 . Consequently, we shall evaluate the threshold of the proposed codes by the density evolution for GL(GF(2), m).

In the binary case, we can predict the asymptotic decoding performance of LDPC codes transmitted over the general MBIOS channels in the large code length limit by density evolution [START_REF] Richardson | Design of capacity-approaching irregular low-density paritycheck codes[END_REF]. Density evolution is possible also for the non-binary LDPC codes [START_REF] Bennatan | On the application of LDPC codes to arbitrary discrete-memoryless channels[END_REF] but computationally intensive and tractable only for the BEC. The analysis for the BEC often helps us to capture the universal properties of LDPC codes.

When the transmission is taken place over the BEC and all-zero codewords are assumed to be sent, the messages, described by probability vectors (p(x)) x∈GF(2 m ) of length 2 m in general, can be reduced to linear subspaces of GF(2) m [START_REF] Rathi | Density Evolution, Threshold and the Stability Condition for non-binary LDPC Codes[END_REF]. To be precise, for each message in the SP decoding algorithm, a subset of {x ∈ GF(2) m | p(x) = 0} forms a linear subspace of GF(2) m , where x is the binary representation of x ∈ GF(2 m ).

For messages in SP decoding, probability vectors P = (P 0 , . . . , P m ) are used for the density evolution and referred to as densities. The i-th entry P i is the probability that a message forms a subspace of dimension i for i = 1, . . . , m. Define two densities P (ℓ) and Q (ℓ) as the densities of messages sent from variable nodes and check nodes at the ℓ-th iteration round, respectively. In [START_REF] Rathi | Non-binary LDPC codes and EXIT like functions[END_REF], Rathi proved that the density that outgoing messages from a variable (resp. check) node of degree 3 with two incoming messages of density P and Q is given by P ⊡ Q (reps. P ⊠ Q). The detail calculation of P ⊡ Q k and P ⊠ Q k are given as follows, respectively,

m i=k k+m-i j=k 2 (i-k)(j-k) i k m -i j -k m j P i Q j , k i=0 k j=k-i 2 (k-i)(k-j) m -i m -k i k -j m m -j P i Q j , where m k = k-1 l=0 2 m -2 l 2 k -2 l is a 2-Gaussian binomial.
Using these two operations of two densities, the density evolution in [START_REF] Rathi | Non-binary LDPC codes and EXIT like functions[END_REF] gives recursive update equations of P (ℓ) and Q (ℓ) for ℓ ≥ 0.

Rathi [START_REF] Rathi | Density Evolution, Threshold and the Stability Condition for non-binary LDPC Codes[END_REF] developed the density evolution for the BEC that tracks probability densities of the dimension of the linear subspaces. For ℓ ≥ 0, the density evolution tracks the probability vectors P (ℓ) and Q (ℓ) which are referred to as densities. The initial messages in Eq. ( 1) can be seen as the intersection of d subspaces of the messages received as the channel outputs.

With overhead ǫ, the decoder has k(1 + ǫ)/C channel outputs transmitted over the channel with capacity C. The number of variable nodes in C 1 is N . It holds that

N = d c d c -2 k m ,
since C 1 is of rate (d c -2)/d c and defined over GF(2 m ). The average number of collected channel outputs per variable node in C 1 is given by

β := (1 + ǫ)m d c -2 d c 1 C . (2)
Denote the number of channel outputs by n := k(1 + ǫ)/C = βN . It follows that the probability R d that a randomly chosen variable node in C 1 has d corresponding channel outputs is given by

R d = n d 1 N d 1 - 1 N n-d . It follows that d≥0 R d x d = 1 N x + 1 - 1 N n = 1 - β(1 -x) n n (n→∞) = e -β(1-x) = d≥0 β d e -β d! x d . (3)
From this, we see the probability that a randomly chosen variable node in C has d corresponding channel outputs in the limit of k → ∞ is

β d e -β d! .
For the binary case, i.e., m = 1, the probability that a variable node has no channel outputs is given by R 0 = e -β . The threshold of the binary (2, d c )regular LDPC codes is given by 1 dc-1 [START_REF] Richardson | Modern Coding Theory[END_REF]. Therefore the asymptotic overhead of the proposed scheme with binary (2, d c )-regular LDPC codes is given as ǫ such that

e -β = 1 d c -1 .
From this and (2), it follows that

ǫ = d c d c -2 ln(d c -1) -1.
The density of the initial messages is given by P (0) as follows,

P (0) = d≥0 β d e -β d! d times E ⊡ • • • ⊡ E,
where E is a density such that the subspace is of dimension m -1 with probability 1. To be precise, E := (E 0 , . . . , E m ) with

E i := 1 if i = m -1 0 if i = m -1.
Since the pre-code is a (2,3)-regular LDPC codes, we have recursive update equations of densities as follows.

Q (ℓ+1) = P (ℓ) ⊠ P (ℓ) , P (ℓ+1) = P (0) ⊡ Q (ℓ+1) .

Since the messages of dimension 0 corresponds to the successful decoding, the asymptotic overhead ǫ * is defined as follows.

ǫ * := inf ǫ∈[0,1] {ǫ ∈ [0, 1] | lim ℓ→∞ P (ℓ) 0 = 1}.
It follows that, in the limit of many information bits k → ∞, with overhead ǫ>ǫ * the reliable transmissions are possible with the proposed C ∞ .

Table I shows the asymptotic overhead ǫ * of the proposed code C ∞ over GF(2 m ) for different m = 1, . . . , 19. Table I also lists the asymptotic overheads with (2,d c )-regular non-binary LDPC pre-code for d c =4, 5 and 6. It can be seen that the best overhead ǫ * = 0.079 is attained at d c = 3 and m = 9 and the fountain code C ∞ exhibit very poor overhead if defined on GF(2 m ) with m = 1, i.e,. the binary field. We will employ m = 8 in the next section for its good asymptotic overhead ǫ * = 0.0809 and friendliness for byte-oriented processors.

A. Monte Carlo Density Evolution

The thresholds for noisy channels (C=0.5 and C=0.1) are computed with Monte Carlo approximation of density evolution (DE). In the case of non-binary codes, it is no longer possible to rely on histograms to get numerical approximation of the thresholds since the densities of the messages are multi-variable functions. The only way to get thresholds is then to turn to Monte Carlo approximation of the densities of messages, that is using samples drawn from the densities. More details on the algorithm used to compute the thresholds can be found in [START_REF] Sassatelli | Nonbinary hybrid LDPC codes[END_REF] In this paper, we have used an extra variance reduction technique in order to get accurate values for the thresholds. A posteriori averaging is used based on 200 threshold computations initialized with different random seeds. The values for C=1.0 indicated in Table II are computed with our Monte Carlo algorithm with variance reduction, and we can see that the thresholds match perfectly the exact thresholds of Table II. We therefore have a good confidence in the threshold values for the other capacities C=0.5 and C=0.1.

V. NUMERICAL RESULTS

In this section, we present demonstrations of C ∞ defined over GF (2 8 ) with moderate number of information bits. Transmission over the BEC and the BIAWGN channels are considered.

Figure 3 shows the histograms of overheads of C ∞ defined over GF (2 8 ). The number of the information bits is set k = 192, 512, 1024, 2048, 8192, 16384 and 32768 from the top to bottom. It can be seen that the overhead is getting concentrated at the asymptotic overhead 0.0809 as increasing k. This was predicted in Section IV.

Figure 4 shows the decoding performance of the proposed fountain code transmitted over the BIAWGN channels with capacity C = 1.0, 0.5 and 0.1. The horizontal axis describes the overhead and the vertical axis describes the block error rate. The performance of best known raptor codes [START_REF] Venkiah | Rate splitting issue for finite length Raptor codes[END_REF], [START_REF] Venkiah | Analyse et optimisation de codes Raptor pour les canaux multicast sans fil[END_REF], [START_REF] Venkiah | Jointly decoded raptor codes: analysis and design for the BIAWGN channel[END_REF] optimized for k = 1024 are drawn for comparison. The proposed codes exhibit the better performance than the best known raptor codes.

VI. CONCLUSION

In this paper, we proposed a new simple fountain coding scheme. The conventional design of fountain codes faced with problems posed in Section I. Despite no optimization of the output degree distribution is needed, the proposed fountain codes outperform the best know optimized raptor codes for moderate code length. Furthermore, the decoding complexity does not depends on channels. Decoding performance of the proposed fountain codes for the BIAWGN channels with capacity C=1.0, 0.5 and 0.1, The information size is k = 1024. The performance of best known raptor codes [START_REF] Venkiah | Rate splitting issue for finite length Raptor codes[END_REF], [START_REF] Venkiah | Analyse et optimisation de codes Raptor pour les canaux multicast sans fil[END_REF], [START_REF] Venkiah | Jointly decoded raptor codes: analysis and design for the BIAWGN channel[END_REF] optimized for k = 1024 are drawn for comparison.

1 .

 1 Each variable node v = 1, . . . , N in C 1 sends the initial message p v ∈ R 2 m to each adjacent check node c. Set the iteration round as ℓ := 0. check to variable output : For each check node c = 1, . . . , M in C 1 , let ∂c be the set of the adjacent variable nodes. It holds that #∂c = 3, since the pre-code C 1 is (2, 3)-regular. Each c has 3 incoming messages p (ℓ) vc for v ∈ ∂c from the 3 adjacent variable nodes. The check node c sends the following message p

Fig. 3 .

 3 Fig. 3. Histograms of the overheads at which the proposed fountain code over GF(2 8 ) successfully recovers k information bits over the channel with C = 1.0. The number of the information bits is set k = 192, 512, 1024, 2048, 8192, 16384 and 32768 from the top to bottom. The horizontal axis describes the overhead ǫ. It can be seen that it is getting concentrated at overhead 0.0809 as predicted in Tab. I at m = 8.

TABLE I ASYMPTOTIC

 I OVERHEAD ǫ * OF THE PROPOSED CODES C∞ WITH A PRE-CODE (2, dc)-REGULAR NON-BINARY LDPC CODES OVER GF(2 m ) TRANSMITTED OVER THE BEC, I.E., C = 1.0. THE BEST OVERHEADS FOR EACH dc ARE WRITTEN IN BOLD FACE TYPE.

TABLE II THRESHOLD

 II VALUES CALCULATED BY MONTE CARLO DENSITY EVOLUTION FOR C = 1.0 (LEFT), C = 0.5 (CENTER), C = 0.1 (RIGHT).

	C = 1.0	C = 0.5
	m dc = 3 dc = 4	m dc = 3 dc = 4
	4 0.2075 0.2665	4 0.2907 0.3582
	6 0.1069 0.1452	6 0.1516 0.1973
	8 0.0809 0.1011	8 0.1150 0.1410
	10 0.0813 0.0856	10 0.1144 0.1214
	12 0.0912 0.0820	12 0.1274 0.1169
	C = 0.1
	m dc = 3 dc = 4
	4 0.4047 0.4816
	6 0.2073 0.2641
	8 0.1376 0.1783
	10 0.1155 0.1408
	12 0.1129 0.1244

PLACE PHOTO HERE

David Declercq was born in June 1971. He graduated his PhD in Statistical Signal Processing 1998. He worked on a new Gaussianity test based on Hermite polynomials properties, and the characterization and the blind identification of non-linear time series. After his PhD, he oriented his researchs towards digital communications, and especially coding theory and iterative decoder design. He started to work on LDPC codes in 1999, both from the code and decoder design aspects.

Since 2003, he focused on studying and developping LDPC codes and decoders in high order Galois fields GF(q), with q ≫ 2. A large part of his research projects are related to non binary LDPC codes. He mainly investigated two directions: (i) the design of GF(q) LDPC codes for short and moderate lengths, and (ii) the simplification of the iterative decoders for GF(q) LDPC codes with complexity/performance tradeoff constraints.

David Declercq published more than 20 papers in major journals (IEEE-Trans. Commun., IEEE-Trans. Inf. Theo., Commun. Letters, EURASIP JWCN), and more than 70 papers in major conferences in ICT. He is currently full professor at the ENSEA in Cergy-Pontoise, France, a graduate school in Electrical Engineering. He is a member of the ETIS laboratory, general secretary of the National GRETSI association, and member of the GdR-ISIS direction team. He is currently the recipient of junior position at the "Institut Universitaire de France".

PLACE PHOTO HERE

Kohichi Sakaniwa received B.E., M.E., and Ph.D. degrees all in electronic engineering from the Tokyo Institute of Technology, Tokyo Japan, in 1972Japan, in , 1974Japan, in and 1977, ,


Jean-Baptiste Fiot

Laurent Risser

Laurent D Cohen

Jurgen Fripp

François-Xavier Vialard

Local vs global descriptors of hippocampus shape evolution for Alzheimer's longitudinal population analysis

Keywords: Brain imaging, population analysis, Alzheimer's disease, geodesic shooting, time-series image data, Kärcher mean

In the context of Alzheimer's disease (AD), state-of-the-art methods separating normal control (NC) from AD patients or CN from progressive MCI (mild cognitive impairment patients converting to AD) achieve decent classification rates. However, they all perform poorly at separating stable MCI (MCI patients not converting to AD) and progressive MCI. Instead of using features extracted from a single temporal point, we address this problem using descriptors of the hippocampus evolutions between two time points. To encode the transformation, we use the framework of large deformations by diffeomorphisms that provides geodesic evolutions. To perform statistics on those local features in a common coordinate system, we introduce an extension of the Kärcher mean algorithm that defines the template modulo rigid registrations, and an initialization criterion that provides a final template leading to better matching with the patients. Finally, as local descriptors transported to this template do not directly perform as well as global descriptors (e.g. volume difference), we propose a novel strategy combining the use of initial momentum from geodesic shooting, extended Kärcher algorithm, density transport and integration on a hippocampus subregion, which is able to outperform global descriptors.

Introduction

Large scale population studies aim to improve the understanding of the causes of diseases, define biomarkers for early diagnosis, and develop preventive treatments. For Alzheimer's disease, several classification strategies have been proposed to separate patients according to their diagnosis. These methods can be split into three categories: voxel-based [START_REF] Klöppel | Automatic classification of MR scans in Alzheimer's disease[END_REF][START_REF] Vemuri | Alzheimer's disease diagnosis in individual subjects using structural MR images: Validation studies[END_REF][START_REF] Lao | Morphological classification of brains via high-dimensional shape transformations and machine learning methods[END_REF][START_REF] Magnin | Support vector machine-based classification of Alzheimer's disease from whole-brain anatomical MRI[END_REF][START_REF] Fan | Compare: classification of morphological patterns using ADaptive regional elements[END_REF][START_REF] Fan | Alzheimer's Disease Neuroimaging Initiative: Spatial patterns of brain atrophy in MCI patients, identified via high-dimensional pattern classification, predict subsequent cognitive decline[END_REF][START_REF] Fan | Structural and functional biomarkers of prodromal Alzheimer's disease: a high-dimensional pattern classification study[END_REF], cortical-thickness-based [START_REF] Klöppel | Automatic classification of MR scans in Alzheimer's disease[END_REF][START_REF] Querbes | Early diagnosis of Alzheimer's disease using cortical thickness: impact of cognitive reserve[END_REF][START_REF] Desikan | Automated MRI measures identify individuals with mild cognitive impairment and Alzheimer's disease[END_REF] and hippocampus-based [START_REF] Chupin | Anatomically constrained region deformation for the automated segmentation of the hippocampus and the amygdala: Method and validation on controls and patients with Alzheimer's disease[END_REF][START_REF] Chupin | Fully automatic hippocampus segmentation and classification in Alzheimer's disease and mild cognitive impairment applied on data from ADNI[END_REF][START_REF] Gerardin | Multidimensional classification of hippocampal shape features discriminates Alzheimer's disease and mild cognitive impairment from normal aging[END_REF] methods. While decent classification rates can be achieved to separate AD from CN or CN from progressive MCI (MCI patients converting to AD), all methods perform poorly at separating stable MCI (MCI patients non converting to AD) and progressive MCI. A recent review comparing these methods can be found in [START_REF] Cuingnet | Automatic classification of patients with Alzheimer's disease from structural MRI: A comparison of ten methods using the ADNI database[END_REF].

In this paper, we investigate the use of longitudinal evolution quantifiers either local or global to separate between stable MCI and progressive MCI. To extract information between two successive time-points, we use a one-to-one deformation mapping the first image onto the second one. Different registration algorithms are available to compute plausible deformations in this context. However, only one, the large deformations via diffeomorphisms (LDD) [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF], provides a Riemannian setting that enables to represent the deformations using tangent vectors: initial velocity fields or equivalently initial momenta. This can be used in practice to retrieve local information and to perform statistics on it as presented in [START_REF] Vaillant | Statistics on diffeomorphisms via tangent space representations[END_REF][START_REF] Wang | Large deformation diffeomorphism and momentum based hippocampal shape discrimination in dementia of the alzheimer type[END_REF]. In this direction, it is worth mentioning paper [START_REF] Singh | Multivariate statistical analysis of deformation momenta relating anatomical shape to neuropsychological measures[END_REF] which shows the correlation between principal modes of deformations and diagnosis. In order to compare this information among the population, we need to define a common coordinate system. This implies (1) the definition of a template and (2) a methodology for the transport of the tangent vector information.

In the literature, point ( 1) is addressed via different methods [START_REF] Ma | A bayesian generative model for surface template estimation[END_REF][START_REF] Fletcher | Principal geodesic analysis for the study of nonlinear statistics of shape[END_REF]. Combining geodesic shooting algorithm presented in [START_REF] Vialard | Diffeomorphic 3D image registration via geodesic shooting using an efficient ADjoint calculation[END_REF], we chose to develop a Kärcher method to average a set of shapes. A first approach has been presented in [START_REF] Vialard | Diffeomorphic 3d image registration via geodesic shooting using an efficient ADjoint calculation[END_REF]. A natural requirement on the Kärcher average is that it could be invariant with respect to rigid transformations of each subject of the population. However, this is not the case in [START_REF] Vialard | Diffeomorphic 3d image registration via geodesic shooting using an efficient ADjoint calculation[END_REF]. One of the contribution of the present paper is to propose a methodology to define such invariant Kärcher averages. We also use a finer strategy to update the deformations at each iteration of the algorithm. Point (2) benefited in each different settings from various contributions that go beyond the standard transport actions. The key-point in our application is that inter-subject variability is much higher than the longitudinal variation so that one expects the statistical results to be strongly influenced by the choice of the transport. To address this issue, parallel transport has been proposed in the LDD setting in [START_REF] Younes | Jacobi fields in groups of diffeomorphisms and applications[END_REF] and it has been applied to longitudinal data discrimination, very similar to our problem, in [START_REF] Qiu | Parallel transport in diffeomorphisms distinguishes the time-dependent pattern of hippocampal surface deformation due to healthy aging and the dementia of the Alzheimer's type[END_REF]. Note that parallel transport preserves the norm of the velocity field and since this norm is not invariant with respect to rescaling, the population variability is still contained in the parallel transport of the tangent information. For other frameworks, such as Log-demons, Schild's ladder approach has been introduced in [START_REF] Lorenzi | Geodesics, parallel transport & one-parameter subgroups for diffeomorphic image registration[END_REF] to extend parallel transport to their Lie-group setting. In any case, we consider the question of transporting tangent information as still open and this motivates us to compare different transport strategies in the classification step.

Section 2 introduces the global pipeline used to perform statistics from local descriptors of hippocampus deformations, with details on geodesic shooting, population template and transport. Section 3 presents the data used and the numerical results. Section 4 concludes the paper.

Methods

Global pipeline

Let us assume we have a population of patients and the binary segmentation of their hippocampus at two different time points. Let us also assume that all patients have the same diagnosis at the first time point, and only a part of them have converted to another diagnosis at the second time point. Our goal is to compare patient evolutions, and classify them with regard to stable diagnosis versus progressive diagnosis.

We use the pipeline summarized in Fig. 1 to compute descriptors of patient evolutions in a common space. First, the evolution descriptors is computed locally for each patient (independently). To be able to compare these descriptors, one needs to transport them in a common space. To do so, a population template is computed, towards which all the local descriptors are transported. Finally, classification is performed to separate progressive from stable patients. Several local descriptors were tested: initial momentum and initial velocity field of geodesic shooting. The use of subregion and integration were also introduced. As for global shape deformation descriptors, volume variation and relative volume variation were tested.

Input images

Local deformation descriptors

Population template

Transported deformation descriptors Classification Fig. 1: Four steps are needed to classify patient evolutions using local descriptors of shape deformations: (1) the local descriptors are computed for each patient independently, (2) a population template is computed, (3) all local shape deformation descriptors are transported towards this template, (4) SVM classification is performed for local descriptors on the whole domain and on a subregion, and local descriptors integrated on the whole domain and on a subregion, and compared with results using global descriptors.

Geodesic shooting

The LDD framework introduces the following minimization problem:

argmin v∈L 2 ([0,1],V ) 1 2 I • φ -1 0,1 -J 2 L 2 + λ 1 0 |v t | 2 dt , ( 1
)
where I is the source image, J the target image, and v t is a time dependent velocity field that belongs to a reproducing kernel Hilbert space of smooth enough vector fields. The deformation Φ is given by the flow of v t . Namely, for any x ∈ Ω the domain of the images,

∂ t φ 0,t (x) = v t (φ 0,t (x)) φ 0,0 (x) = x .
The first term in formula ( 1) is a similarity term controlling the matching quality whereas the second one is a smoothing term controlling the deformation regularity. The minimization problem (1) can be reformulated using a shooting formulation on an initial condition which is the so-called initial momentum denoted P (0), a scalar field in the image space, as follows:

argmin

P (0) 1 2 I • φ -1 0,1 -J 2 L 2 + λ ∇I(0)P (0), K ⋆ ∇I(0)P (0) L 2 (2)
subject to the geodesic shooting equations:

     ∂ t I + ∇I • v = 0 ∂ t P + ∇ • (P v) = 0 v t + K ⋆ ∇IP = 0 , (3)
where I is the source image, v t the velocity field, P the momentum, K the translation-invariant kernel of the reproducing kernel Hilbert space, and ⋆ the convolution operator. In order to solve the new optimization problem (2), we use the methodology described in [START_REF] Vialard | Diffeomorphic 3D image registration via geodesic shooting using an efficient ADjoint calculation[END_REF]. Note that the choice of the kernel matters for retrieving plausible deformations and we refer to [START_REF] Risser | Simultaneous multiscale registration using large deformation diffeomorphic metric mapping[END_REF] for an extensive discussion on the parameter choices.

For each patient, a two-step process was performed to encode the deformations of the hippocampus shape evolution, as described in Fig. 2. First the hippocampus from the second time point was rigidly registered back to the first time point. Second, the geodesic shooting was performed from the first time point towards the registered second time point. Inital momenta and initial velocity fields from different patients are local descriptors that were used to compare hippocampus evolutions.

Population template

As mentioned in section 2.1, local descriptors of hippocampus evolutions need to be transported in a common space prior to any statistical analysis. We introduce two modifications of the standard Kärcher mean algorithm [START_REF] Vialard | Diffeomorphic 3d image registration via geodesic shooting using an efficient ADjoint calculation[END_REF]: (1) the population template is computed up to rigid transformations and (2) the template is regenerated from a reference image at every iteration. We propose the following optimization problem:

K = argmin K,R1,...,Rn n i=1 d(K, I i • R -1 i ) 2 , ( 4
)
Step 1. Rigid registration

Step 2. Shooting Input:

source image: J (follow-up, t = t0 + 12 months) -target image: I (screening, t = t0) Output:

rigid transformation R Input:

source image: I target image: R(J) Output:

initial momentum P 0 Fig. 2: For each patient, the initial momentum encoding the hippocampus evolution is computed in a two-step process. First, the second time point (t = 12 months) is rigidly registered to the first time point (t = 0). Second, the geodesic shooting is computed from the first time point to the previously rigidly registered second time point.

where {I 1 , . . . , I n } are the population images, {R 1 , . . . , R n } are the rigid transformations registering the I i towards K.

The proposed associated algorithm initializes the deformation field generating the template from the reference as the identity: u 0 = Id. Each Kärcher iteration is then composed of four steps described in Fig. 3.

The first modification is introduced because current implementations of LDD algorithms are not invariant by the action of the group of orthogonal linear transformations, so the resulting template reflects the orientation variability of the population. Numerically speaking, such a modification can potentially help avoiding convergence issues of the geodesic shootings (if the template estimate and the image are too far with respect to the LDD metric), and reduce computation time. The second modification of the standard Kärcher algorithm is purely motivated by numerical considerations. The standard formulation smoothes the template estimate at each iteration (as K t+1 is computed by shooting from K t ). In our algorithm, the template is regenerated from the reference image at each step, so that the template estimate keeps sharp boundaries.

Tangent information and associated transport

The local descriptors computed for each patient as explained in section 2.2 need to be transported in a common coordinate space: the space of the Kärcher average defined in section 2.3. We chose transport rules that only depend on the final deformation (it does not include parallel transport which depends on the chosen path). A two-step process was then used to transport local descriptors Fig. 3: Each Kärcher iteration is composed of four steps: (1) the images are rigidly aligned towards the current Kärcher mean estimate, (2) geodesic shootings from the current Kärcher estimate K t towards all the registered images are computed (3) geodesic shooting from K t using P 0 = 1 n i P 0 i generates a deformation field u tmp , and (4) the composed deformation field u t+1 = u tmp • u t is used to compute the updated estimate from the reference image. of hippocampus evolutions to the template space (Fig. 4). First, the screening hippocampus was registered towards the template rigidly [START_REF] Ourselin | Reconstructing a 3D structure from serial histological sections[END_REF] then non-rigidly [START_REF] Modat | Fast free-form deformation using graphics processing units[END_REF]. The resulting deformation is denoted by φ. Second, this transformation was used to transport the local descriptors of hippocampus deformations towards the template. The transport itself depends on the nature of the quantity transported. For instance, we call image transport the standard transformation of an image by the action of a diffeomoprhism φ:

I transported (x) = I local • φ -1 (x) . (5)
From the mathematical point of view, the momentum is an adjoint variable to the image. As a consequence, it is transported by the adjoint action of the group which reduces to the standard transport for a density n, defined by:

n transported (x) = det(Jac φ -1 (x)) n local • φ -1 (x) . ( 6
)
On the theoretical side, this transport preserves geodesics, due to the Ad invariance of momentum maps: we refer the reader to [START_REF] Holm | Geometric Mechanics and Symmetry: From Finite to Infinite Dimensions[END_REF] for that property. Note that this action preserves the global integration of the density by a simple change of variable. Last, we present the transport via the standard conjugation of a velocity field defined by:

V transported (x) = dφ(x) • V local • φ -1 (x) . (7)
All those transport methods were tested in the classification step. We did not include parallel transport in this study since no public implementation is available and its implementation is rather involved, especially in the case of images. 3 Material and Results

Data

A dataset of 206 hippocampus binary segmentations from 103 patients enrolled in ADNI1 [START_REF] Mueller | The Alzheimer's Disease Neuroimaging Initiative[END_REF] has been used to estimate the efficiency of local and global descriptors of hippocampus evolution with regard to disease progression. For each patient, 'screening' and 'month 12' were the two time points selected. All patients were MCI at the screening point, 19 became AD by month 12, and the remaining 84 stayed MCI.

Experiments

First, all screening images were resampled to a common isotropic voxel size 1.0 × 1.0 × 1.0 mm, similar to their original size. Rigid transformations aligning the month 12 hippocampus towards the screening ones were computed using [START_REF] Ourselin | Reconstructing a 3D structure from serial histological sections[END_REF]. The geodesic shootings [START_REF] Vialard | Diffeomorphic 3D image registration via geodesic shooting using an efficient ADjoint calculation[END_REF] were performed2 using a sum of three kernels (sizes 1, 3 and 6 mm, with respective weights 2, 1 and 1), and 200 gradient descent iterations. All the 103 source images deformed by the shooting had a final Dice score (i.e. overlap score) with their target images over 0.85.

To compute the template, a subset of 20 images was used. This subset and the initialization was based on a shape volume criterion. Four Kärcher iterations were performed, with respectively 200, 150, 150 and 100 gradient descent iterations in the geodesic shootings. To compute the transformations from the screening hippocampi towards the template, rigid [START_REF] Ourselin | Reconstructing a 3D structure from serial histological sections[END_REF] then non-rigid [START_REF] Modat | Fast free-form deformation using graphics processing units[END_REF] registration algorithms were applied with their default parameters.

To classify from local descriptors, a mask computed by dilating the template was used. To compute classification on subregions, each hippocampus (left and right) from the template was dilated. The bounding box was cut equally in thirds along the longest axis, and intersections were used as masks.

Using a leave-10%-out scheme, training and test sets were created. With training features equally distributed among classes, SVM classifiers were computed (the Matlab functions from the Bioinformatics Toolbox were used). All the patients were then classified. The Gaussian kernel was used and 20 kernel widths were tested. This procedure was repeated 50 times and classification accuracy averaged. From the numbers of true/false positives/negatives (TP, FP, TN, FN), four indicators were used to measure classification accuracy: specificity Spec = T N T N +F P , sensitivity Sens = T P T P +F N , negative predictive value N P V = T N T N +F N , and positive predictive value P P V = T P T P +F P .

Results

Using the modified Kärcher mean algorithm and the criterion mentioned above, the average Dice score between the 103 registered patients and the template was 0.87 ± 0.02, whereas it was only 0.44 ± 0.11 when matching to a template computed using a criterion based on the distance to the L2 mean.

Regarding descriptors of hippocampus evolutions, the local descriptors did not perform as well as global descriptors, when used directly as input features (Fig. 5). However, when integrated on the whole domain, the performances were similar. When integrated on some subregion, they can outperform the global descriptors. Detailed results are displayed in Fig. 5, and Table 1 displays the four unbiased indicators when the sum of specificity and sensitivity is maximized.

Discussion and conclusion

We have studied the use of global, semi-local and local descriptors of hippocampus evolutions to predict AD conversion for MCI patients, using a dataset of binary segmentations provided by ADNI. This study focuses on shape evolutions between two time points, whereas (to the best of our knowledge), studies in this application field usually extract features from a single time point and perform diagnosis classification.

The proposed extension of the Kärcher mean algorithm with a subpopulation and initialization criterion based on shape volume improved the matching quality to the template (average Dice of 0.87 ± 0.02 instead of 0.44 ± 0.11) without the need of modifying the default registration parameters. In our experiments, the local descriptors did not perform as well as global descriptors such as volume difference when they were directly used as input features of the SVM classifiers. However, when integrated over the whole domain, classification performances were similar. When integrated on a subregion, they could even outperform the global descriptors. The method we propose combines (1) the use of initial momentum of geodesic shooting, (2) an extended version of the Kärcher mean algorithm, (3) the use of density transport and (4) the integration on a subregion. On our dataset, this method was the only one able to outperform the global descriptors. It should be noted that in our study the definition of the subregion was sub-optimal and used as a proof-of-concept. The most promising perspectives are (1) developing a strategy to define subregions maximizing the classification results and (2) adding more time-points to the study using the geodesic regression method introduced in [START_REF] Niethammer | Geodesic regression for image timeseries[END_REF] or cubic spline interpolation in [START_REF] Trouvé | Shape splines and stochastic shape evolutions: A second order point of view[END_REF].

Fig. 4 :

 4 Fig.4: Local descriptors of hippocampus evolutions are transported to the template in a two-step process. First the deformation field from the patient space to the population template. Second, this deformation field is used to transport the local descriptors.

Fig. 5 :

 5 Fig.5: Classification performance (depending on the SVM Gaussian kernel width σ) for global descriptors (5a, 5b), local (5c, 5f, 5i), local integrated on the whole image (5d, 5g, 5j) and local integrated on a subregion (5e, 5h, 5k). Higher for Spec + Sens (in cyan blue) is better.

Table 1 :

 1 Performance indicators for various descriptors of the hippocampus evolutions. These indicators are computed using a SVM classifier, with a Gaussian kernel. Kernel width is such that the sum of specificity and sensitivity is maximized. The proposed method is the only one with Spec + Sens outperforming the same sum for the volume difference global descriptor. Initial momentum, image transport 1.18 0.63 0.55 0.86 0.25 Initial momentum, density transport 1.27 0.62 0.65 0.89 0.28 Initial momentum, image transport 1.10 0.68 0.42 0.84 0.23 Initial momentum, density transport 1.17 0.68 0.49 0.85 0.26 Initial velocity field 0.98 0.38 0.60 0.81 0.18

	Global / Local	Deformation descriptor	Spec+ Spec Sens NPV PPV Sens
	Global	Volume difference Relative volume difference	1.19 0.78 0.41 0.85 0.30 1.08 0.85 0.23 0.83 0.25
	Local integrated	Initial momentum, image transport 1.10 0.37 0.73 0.86 0.21
	on the whole	Initial momentum, density transport 1.15 0.96 0.19 0.84 0.53
	domain	Initial velocity field	1.07 0.46 0.61 0.84 0.20
	Local integrated
	on a subregion	Initial velocity field	0.92 0.79 0.13 0.80 0.11
		Initial momentum, image transport 1.01 0.96 0.05 0.82 0.27
	Local	Initial momentum, density transport 1.01 0.95 0.06 0.82 0.21
		Initial velocity field	0.92 0.77 0.15 0.80 0.13
	Local restricted
	to a subregion

http://www.loni.ucla.edu/ADNI

http://sourceforge.net/projects/utilzreg/


Philippe Lambert

Sylviane R Schwer

Nicolas Boffo

A new model of Time Expressions Detection and Annotation in Vietnamese : The hôm case

Keywords: Natural Language Processing, automatic translation (Vietnamese-French), Vietnamese calendrical expressions, temporal annotation and reasoning, Viet4Nooj module

We describ e our approach to build a new model of Time Expressions Recognition in Vietna mese, in cluding a temporal reasoning tool in order to provid e an automatic sup ervised lea rning system based on a socio-cultural linguisti c approach. The Vietna mese morpheme hôm is taken for examplification of th e impediments to be solved .

I. INTRODUCTION

The importance of modeling temporal information is increasingly apparent in natural language applications, such as information extraction, question answering, translation or second language acquisition. Taking into account temporality expressed in texts appears as fundamental not only in a perspective of global processing of documents but also in the comprehension of the structure of a document. A lot of works have put the emphasis on the importance of the temporal expressions as modes of discursive organization. Among temporal expressions, relative calendrical expressions such as Sunday, today, yesterday morning, the day after tomorrow, three days before can be viewed as a part of Named Entity Resolution -they refer to calendrical periods -strongly depending on contextual or textual cues. In learning and translating perspective, the usual practice of encoding the calendrical periods designating by temporal expressions in terms of dates is far to be sufficient. Our Claim is that different expressions referring to the same period design different perspectives of the writers, which have to be kept into account during the comprehension and the translation process. For instance, Fr. Pâques 2012 refers to the same date as le 8 avril 2012 but the former refers to a festival day with all its religious characters, the latter to a simple Gregorian date. The same remark can be done in the Vietnamese context with the celebration doctors day (Ngày Th1y thu2c VN, February 27) or the Liberation Day of Sài Gòn (Ngày Gi3i phóng Mi4n Nam, April 30) . Some calendrical expression can have different meanings, depending on co-text or pragmatic situation. But none are able to distinguish the different meanings of Fr. hier/aujourd'hui/demain, that can either refer to the day-unit [START_REF] Puste Jovs Ky | Ti meBank 1 . 2 Doc u menta tion[END_REF] or to the partition of a wider span of time [START_REF] Do | Extracti on de c or pus para llè le pour la traducti on au toma tique de puis e t vers une langue peu dotée[END_REF].

( The hôm case.

language TimeML in [START_REF] Puste Jovs Ky | Ti meBank 1 . 2 Doc u menta tion[END_REF] next Tuesday is encapsulated inside the normalized form "two weeks from" and provides the Gregorian date associated as shown in (3).

(3) <T IMEX2 VAL ="1999-08-03">two weeks from <T IMEX2 VAL="1999-07-20">next Tuesday</ TIMEX2 ></TIMEX2> Actually, the automatic translation systems are based on probabilitic and statistic functions. For instance, they have difficulty in differentiating between hôm and ngày, as shown in (4). Hôm nay cannot be associated in this context to Aujourd'hui, which is used as in ex. [START_REF] Do | Extracti on de c or pus para llè le pour la traducti on au toma tique de puis e t vers une langue peu dotée[END_REF]. One has to choose between ngày nay or hiAn t7i. In [START_REF] Do | Extracti on de c or pus para llè le pour la traducti on au toma tique de puis e t vers une langue peu dotée[END_REF] is mentioned the difficulty in dealing correctly with temporal expressions inside probabilistic and statistic translation systems. In Natural Language Processing community, the first implementation of some Vietnamese temporal expressions has been made inside the NooJ language environment [START_REF] Ilber Zte In | N ooJ : a li nguis tic annota t ion s ys te m for c orpus pr ocess ing[END_REF]. This language is based on final automata transducers, or other said, on formal regular grammars. The Viet4NooJ (V4N) module refers to the linguistic resources created specifically for data treatments in Vietnamese [START_REF] La | A Vie tna mese modu le for NooJ : M ode lizati on, r ea liza tion and perspec tives[END_REF][START_REF] La Mbert | V IET4Nooj A V ie tna mese mo dule for Nooj[END_REF]. A first attempt to add a temporal component has been made for dates' extraction [START_REF] Inh | Au toma tic Process ing of Te mpora lity for VIET4Nooj[END_REF]. We will get advantage of the very beginning of this task to create a multidiscipline group around temporality in Vietnamese, in order to build a complete temporal system, based on deep linguistic and socio-cultural analysis. In the following, we outline the framework in which we wish to contribute in one hand to the resolution of these impediments and, in the other hand, to enrich the annotation in order to include a temporal reasoning system. In the next section, we give a general survey of our multilevel annotations description. We develop the interest -particularly for a temporal reasoning task -of choosing the S[et]languages framework [START_REF] Schwer | Représentation mathématique du temps : après Reichenbach[END_REF] for the annotation of temporal relations provided by deictic calendrical expressions. Then we present the implementation part.

II. ABOUT ANNOTATIONS

Based on the presupposition that two expressions that refer to the same object differ on what they mean about the situation of the discourse (either aspectual, culturally, intentionally, and so on), we propose a multilevel annotation description in two parts: one for the reference and one for the situation of discourse (in large). We aim at providing a general framework of annotation taking into account the two parts of the description and the possibility of interconnection between other annotating frameworks such as TimeML and the temporal reasoning tool SLS developed by Irène Durand in labri, based on the Sylviane S. Schwer's S-languages framework [START_REF] Durand | A Tool for Reasoning about Qualitative Temporal Information: the Theory of Slanguages with a Lisp Implementation[END_REF]. The first step, in development, concerns the calendrical unit "day", which is the basic and more complex one, due to the fact that in one hand, this is the canonical measure of time, and or the other hand, the day unit is based on the most fundamental light/night alternative which still determines (in some extend) activities and rests of human beings.

The temporal reference part: This part can be viewed as the classification of all recognizing calendrical expressions containing these three words according to their denotation. For instance, the three Viet. expressions hôm nay, bBa hôm nay, ngày hôm nay denote the same deictic relationship: they relate some processes inside a unitday-period that contains the moment of speech (coded S). bBa nay and ngày nay denote another deictic relationship: they relate some processes inside a period (to be further specified) that contains the moment of speech. The Fr. Aujourd'hui subsumes the five preceding Vietnamese expressions, and could be described in term of a period. In order to deal also with temporal reasoning, these relationships are encoded in terms of S-languages, which can be viewed as an algebraic formal language that generalizes the reichenbachian way of coding temporal relations expressed by tense [START_REF] Ichenbach | Ele ments of Symbolic Logic[END_REF]. Let us go on with calendrical deictic expressions, that is relative to the moment of speech.

Each calendar unit is associated to a letter (ex. d for day, y for year). An occurrence of a letter depicts a bound (the | in Figure 1) between two successive durations of the unit.

- -----|-----|---E'---|---1--E-|------|------|---E"----| Figure 1 For instance, today is described as +d[S⊗?]d+. I n t h i s e x p r e s s i o n , "?" is a letter that signals something is said to occur inside the duration defined by dd and "[S⊗?]" says that the expression gives no information that allow to order S and ?. The exact temporal relation between S and E is given either by the cotext or by the situation of discourse.

+d and d+ signify that d is iterable both backward and forward.

+d[S⊗E]d+ encodes the fact that event E occurred/occurs/will occur today. Yesterday is described as +d?dSd+, tomorrow is described as +dSd?d+, ten days ago is described as +d?d 10 Sd+, this year is described as +y[S⊗?]y+, ten years ago is described as +y?y 10 Sy+.

Let the sentence (5) be an illustration of Figure1, (5) Yesterday I saw Eve (E'), today I will see Jean (E) and in three days Zoé (E"). its S-expression is the result of the Join operation +dE'dSd+, +d[S⊗E]d+, SE1 ,+dSd 3 E"d+, that is: +dE'dSEd 3 E"d, which is depicted Figure 1, where | stands for d. SLS [START_REF] Schwer | Représentation mathématique du temps : après Reichenbach[END_REF] is a software that computes the Join of any set of S-expressions. An interface between Nooj and SLS is intended.

Based upon the S-language annotation, it is trivial to see that +d?dSd+ and +y?dSy+ can be subsumed by the schema +u?uSu+, that describes that something is said on the unit preceding the deictic anchoring unit encoded by u. This allows a hierarchical structure of expressions, not based on the kind of unit, but only on any temporal partitioning and the proximity to the deictic anchor. For instance, u can encoded any saillant type of events such that a sunami occurrence: in that case +u?uSu+ codes a process that occurred between the the sumami that preceded the last one and the last one. This is also the good way to encode the Fr deictic hier/aujourd'hui/demain, that is +u?uSu+, +uS⊗? ? ?

?u+ , +uSu?u+ with u=d as in [START_REF] Puste Jovs Ky | Ti meBank 1 . 2 Doc u menta tion[END_REF] and u= change of period, to be determined as political phase of a state in [START_REF] Do | Extracti on de c or pus para llè le pour la traducti on au toma tique de puis e t vers une langue peu dotée[END_REF].

The contextual part: this is the sociolinguistic analysis part. The discrimination will be made between each expression of the same class. This requires a deep understanding of the language, in all its aspects. Vietnamese belongs to the group of isolating, tonal a n d m o n o s yl l a b i c languages. Its word forms are invariant, there is no verbal tense. Hence, all grammatical relations are expressed by word order and/or function words and/or tonal and prosody device. Temporality is strongly affected by these characteristics. All linguistic domains have to be questioned.

• Lexically: Northern Vietnamese uses ngày while Southerners prefer bBa, • Syntactically: the position of the temporal adverbial in a question/answer indicates the future or the past : Khi nào anh v4 (future) vs Anh v4 khi nào (past), • Semantically: tone has a semantic value and can be used to describe proximity/remotness: hôm kia/kià/… [START_REF] Nguy | 5n Questions de linguistique vietnamienne : les classificateurs et les déictiques[END_REF],

• Pragmatically: cultural and social aspects linked to calendrical systems and time perception : âm lCch / dDEng . Beyond these classic ones, the following aspects can't be neglected. Vietnamese NLPA studies such as [START_REF] Mac | Modeling the prosody of Vietnamese attitudes for expressive speech synthesis[END_REF][START_REF] Phuong | Finitestate description of Vietnamese reduplication IN[END_REF] testify of their importance.

• Tonaly: specific tonal transformation in regional dialect describing a temporal proximity: in Southern dialect, "hFm" is equivalent to "hôm y".

• Prosody: (i) the use of reduplication can be used for vagueness: chi4u chi4u "epoq epoq"for "n ear this afternoon" or for iteration ngày ngày "jour jour" for "each day"; (ii) the preference for the bisyllabilism in numbering, which obliges to add mùng before numbers 1 to 10. The Vietnamese study of temporality is generally centered on phenomena related to the predicate, particularly tense and aspect and adverbs. Temporal indicators that lie outside this category, most of which involve calendrical terms, have been quite ignored [START_REF] Durand | A Tool for Reasoning about Qualitative Temporal Information: the Theory of Slanguages with a Lisp Implementation[END_REF], [START_REF] Nguy | 5n Questions de linguistique vietnamienne : les classificateurs et les déictiques[END_REF], or [START_REF] Mac | Modeling the prosody of Vietnamese attitudes for expressive speech synthesis[END_REF][START_REF] Phuong | Finitestate description of Vietnamese reduplication IN[END_REF]. For modeling the temporal elements in Vietnamese we proceed as a second language learner who want to understand how to construct a sentence with temporal marks. So we look in many grammatical [START_REF] Nguy | 6n, Nghiên cu ngB pháp ti8ng ViAt[END_REF][START_REF]Boàn Thi9n ThuCt, Nguy6n Khánh Hà, PhDm NhE QuFnh, A concise vietnamese grammar for non-native speakers[END_REF] and pedagogical books [START_REF] Di9p | NgA pháp ti1ng Vi9t [Grammaire vietnamienne[END_REF][START_REF] Lê | Le parler vietnamien[END_REF][START_REF]Ngy6n Vi9t HE4ng, Ti1ng vi9t nâng cao dành cho ngE4i ngoài (quyn 1)[END_REF], dictionaries [START_REF]T in Ti1ng Vi9t, Vietlex, tr ung tâm t in hc, Nhà xut bn Bà Nng[END_REF][START_REF]T in Vi9t Pháp, Vi9n Khoa hc xã h7i Vi9t Nam, Nhà xut bn vn hóa Sài Gòn[END_REF] and lexical [START_REF]Ti1ng vi9t cho ngE4i nE2c ngoài, Nhà Xut B n Vn Hóa[END_REF] of Vietnamese for non-native speakers. Despite their useful analysis, the complements of tense and temporal frames are not completely mentioned and well studied. These works often focus on communication skills without providing an "immediate constituent analysis", which would be helpful for our analysis. For instance, in [START_REF] La Mbert | V IET4Nooj A V ie tna mese mo dule for Nooj[END_REF], "today" is translated by five expressions:

(6) bBa nay, hôm nay, ngày nay, bBa hôm nay, ngày hôm nay. But there is no explanation about the meaning and the syntactical combinatory possibilities of hôm, ngày, and bBa. Actually the automatic translation systems have difficulty in differentiating the constructions with «hôm » and « ngày », expecially when the co(n)text can't be found in their table-phrases, as shown in example [START_REF] La | A Vie tna mese modu le for NooJ : M ode lizati on, r ea liza tion and perspec tives[END_REF]. Some explanations can be found in [START_REF] Kh | T in Pháp Vi9t, Nhá xut bn Thành Ph H Chí Minh[END_REF], where is said that bBa is used in Saigon, "as one of a series, or time when, primarily in present or future", hôm "as time when, primarily in present or past", and ngày "as one of a series, or time when, in future", Today is translated by hôm nay in Hanoi and by bBa nay in Saigon. Nowadays, it seems that hôm nay is used in both area to mean the calendrical present day, bBa hôm nay only in Saigon and ngày hôm nay in Hanoi. BBa hôm nay and ngày hôm nay are used for today/nowadays. Following [START_REF] Th Ompson | A Viet na mese Gr a mmar[END_REF] and [START_REF] Schwe R | Au jour d 'a ujourd[END_REF], we can argue that the trisyllabic expressions maintain the calendrical today inside its calendrical serie, that is, in constrast with yersterday and tomorrow. These explanations provide a reason to have the five different expressions given in [START_REF] Inh | Au toma tic Process ing of Te mpora lity for VIET4Nooj[END_REF] for translating "today" with a dayunit term and the way to differenciate them? III. A FIRST GRAMMAR FOR HÔM AND ITS IMPLEMENTATION This paper constitutes a first work about Vietnamese temporality detection and formalization. It consists not only in providing results of textual treatment but also to test the robustness of detecting temporality structures in texts. Our goal is to provide a set of annoted grammars for each calendrical noun. We have begun with day-unit terms and hôm was choosen as a test, because it is in opposition with the two others. A 1516 articles c o r p o r a w a s b u i l d from the «Th1 Gi2i » category of the Vietnamese website VnExpress.net. We detected then 616 occurrences of «hôm» by using the NooJ's frequency functionality. NooJ 3 is a language environment developed by Max Silberztein [START_REF] Ilber Zte In | N ooJ : a li nguis tic annota t ion s ys te m for c orpus pr ocess ing[END_REF] with a wide open community having developped a set of linguistic tools for more than twenty languages. As we mentioned it above, the Viet4NooJ (V4N) project refers to the linguistic resources created specifically for data treatment in Vietnamese [START_REF] La | A Vie tna mese modu le for NooJ : M ode lizati on, r ea liza tion and perspec tives[END_REF], [START_REF] La Mbert | V IET4Nooj A V ie tna mese mo dule for Nooj[END_REF], [START_REF] Inh | Au toma tic Process ing of Te mpora lity for VIET4Nooj[END_REF]. A few reasons explain why NooJ has been chosen for our work relative to its specificities: Flexibility: NooJ's functionalities are is final end-user oriented plateform. The interface was elaborated to facilitate the creation of linguistic resources for non NLP specialist. Granularity: NooJ uses three particularly efficient computational devices to process texts: (i) Finite-State Automata (FSA) to locate specific expression in corpora, (ii) Recursive Transition Networks (RTN), a grammar coomposed by a set of sub-graphs and (iii) Enhanced Recursive Transition Networks (ERTNs), a set of RTN performing textual operations through the use of variables, producing specific outputs. Interoperability: ERTNs' results (i.e. output / tags) can be re-used in external software. Thus, NooJ would be able to play the role of intermediary phase between heterogeneous copora and computational calculation of time relations with [START_REF] Durand | A Tool for Reasoning about Qualitative Temporal Information: the Theory of Slanguages with a Lisp Implementation[END_REF] as final results. Multimodal tagging process: the system is able to produce specific tags (i.e. created by the user) or xml tags as well. This permits to produce TimeML tags as previously mentioned. Studies about temporal expressions recognition can be schematically described under two main logics : the first one concerns a pattern matching rules approach and the second one, through supervised learning methods [START_REF] Sana Mpudi | Te mpor a l Reason ing in Na tur al Language Pr oc ess ing: A Surve[END_REF]. The semantic based models like Timex2, Chronos or TimeML are the most common systems used by the researchers community working on temporal expression. These models annote temporal expression with two steps: (i) marking a temporal expression in a document, and (ii) identifying the time value that the expression designates. Studies are, until now, focussed on English and european languages such French, German, Italian, slovenian and serbian. Since 2005, studies on Chinese and Korean begun to be developping essentialy to normalize temporal expression and provide tests for both translation and alignement [START_REF] Lecuit | Te mpor a l E xpres s ions : Comparis ons in a M ulti lingua l Cor pus[END_REF][START_REF] Wu | Nor ma lizin g Chinese Te mpora l Expres s ions w ith Mu ltilabe l C las s ifica tion[END_REF]. Concerning more specifically hôm, FSA and RTN were organized to produce a specific ERTNs aiming at detecting the maximal length expression with hôm as core element (Fig. 2). Morphology of Vietnamese language is monosyllabic, every morphem being composed by only one vowel or vocalic component. Graph and subgraph of our ERTNs follow a semantic logic instead of a lexematic one. For example, the expression "hôm th ba" [hôm : day, th : ordinal marker, ba : numeral "3" -> th ba : tuesday] with 3 morphological elements will be tagged with only two hierarchical levels : <hôm,N0><th3 ba,N+1>. This ERTNs is constituted by a set of ten subgraphs organized around the "hôm" nucleus and its regional variant "hFm". The ERTNs is actually able to recognize temporal expression with hôm with a hierachical positioning levels of N-3 -> /.../ N0 (hôm) /->/.../ N+3 The Table 1 gives an example of a relevant expression extracted by V4N: IV. RESULTS AND DISCUSSION Application of the "hôm" ERTNs over the corpora provides good results since the return rate approximates 99 %. However, this result has to be relativized : the test was only conducted on an electronic press sample. We assume that results would have been different with a literary sample. This new test will be s4n explored to pr4ve the robustness of our grammar and to enhance it by adding new lexems to the different levels of the grammar. Beside this, our work can be considered as an opening to some new perspectives on linguistics and NLP for Vietnamese since, as far as the state of art has shown it, till now, almost all of the temporality studies has been focused on temporal lexicography level rather than on detection of temporal framework. Our approach is also bearing a few secondary emerging topics as enhancement of praxematic process, psychocognitive studies, didactic process for non native learners to better understand language structuration and reasoning logic, cross-cultural information retrieval, etc.

V. REFERENCES

( 4 )

 4 Google translation t ranslate the Fr. source "Aujourd'hui nous sommes dans une situation économique difficile." int o th e Viet. Target "Hôm nay chúng ta 5ang 6 trong tìn h tr7ng kinh t8 khó kh9n".

Figure 1 :

 1 Figure 1 : The ER TNs for Hôm

  Very early in the morning on Tuesday

Table 1 :

 1 Hierarchical levels of Hôm : N-3/ 0/N+1

1) Hier la pluie, aujourd'hui la

	A new model of Time Expressions Detection
	and Annotation in Vietnamese :
	pluie et
	demain encore la pluie. (oral sentence,
	recently eared)
	(2) Ainsi, hier une république,
	aujourd'hui une monarchie, et demain
	encore une république.
	CHATEAUBRIAND, Essai sur les
	Révolutions, t. 2, 1797, p. 95.
	Actual time expressions detection and
	annotation systems for different
	languages aim at computing the period
	refered by the expressions and provide
	normalizing expressions as annotations
	[see 1, 2, 3 among others]. For instance,
	using the ISO standard temporal markup

This c o-textua l c onstrain t is given by the fu ture tense .


L Saab

P Souères

N Mansard

J Y Fourquet

Generation of Human-like Motion on Anthropomorphic Systems using Inverse Dynamics

Keywords: humanoid robot, whole-body movement, dynamics, stack of tasks, torque control

Introduction

This work deals with the generation of human-like whole-body movements on anthropomorphic systems. We propose a general framework to generate robot movements from the definition of ordered stack of tasks and a global resolution scheme that enables to consider different kinds of constraints [START_REF] Saab | Generic Dynamic Motion Generation with Multiple Unilateral Constraints[END_REF]. We compare qualitatively the robot movements generated from this software with similar recorded human movements. We start with a direct global comparison of body movements. Then we analyze the magnitude of the reconstructed human torques and compare with the simulated robot torques during the motion.

Humanoid motion generation and human movement analysis

The motion generation method is based on the taskfunction formalism [START_REF] Samson | Robot Control: the Task Function Approach[END_REF]. Complex movements are defined by a sequence of hierarchized stack of tasks [START_REF] Mansard | A Versatile Generalized Inverted Kinematics Implementation for Collaborative Working Humanoid Robots: The Stack of Tasks[END_REF]. Equality tasks are used to drive a task-function to a reference value, whereas inequality constraints enable to maintain functions of state variables above or below a given threshold. The degrees of freedom (DOF) of the robot are then simultaneously constrained to guarantee the execution of the highest priority tasks while satisfying at best the lower priority objectives. All along the movement the robot balance is ensured by specifc constraints.

At each new phase of the movement the stack of task is restructured by introducing or cancelling tasks and by modifying their priority.

A global solver involving a cascade of Quadratic Programming problems [START_REF] Fletcher | A general quadratic programming algorithm[END_REF][START_REF] Escande | F a s t Resolution of Hierarchized Inverse Kinematics with Inequality constraints[END_REF] is used compute the controller. A detailed presentation of this motion generation method is proposed in [START_REF] Saab | Generic Dynamic Motion Generation with Multiple Unilateral Constraints[END_REF]. Note that this framework enables to take into account many possible rigid contacts between the robot and the environment during the movement. Furthermore, this method is generic and can be applied both to the kinematics and to the dynamics of the robot. In this work, we considered an extended model of the humanoid robot HRP-2.

The analysis of the human movement was done at the Yoshihiko Nakamura Lab. (YNL) of the University of Tokyo in collaboration with the LAAS-CNRS. We used a Motion Analysis motion capture device, EMG and force plates to record the movement of human subjects executing the same sequence of reaching movements. Based on these data, the software sDIMS, developed at the YNL, which includes a complete musculoskeletal (MS) model of the human body, was used to reconstruct the torques produced by the subject [START_REF] Nakamura | Somatosensory computation for man-machine interface from motion-capture data and musculoskeletal human model[END_REF]. The main difficulty is to design the ordered stack of tasks to execute the expected sequence of reaching tasks. Figure 2 shows the chronology of the stack of tasks that was considered. The highest priority task ensures the robot balance by maintaining its center of mass (CoM) within the support polygon. Two inequality constraints are used to limit the variation of the left foot and the left hand position during the whole movement while their orientation are free. Three successive equality constraints are used to control the movement of the right hand toward each target. Equality tasks are also used to drive the gaze toward the current target. To avoid auto-collisions and to increase the resemblance to human movement, we also used equality tasks to drive the robot movement through two intermediary waypoints, IP1 and IP2, before reaching the targets. These points were s e l e c t e d a s i n t e r m e d i a r y objectives because they appear to play a key role in the organization of the human motion.

Results

The human m o t i o n a n d t h e s i m u l a t e d r o b o t movements were carefully analyzed. We first observed the global shape of movements and compared them. For this study we developed a specific robot model including 37 DOFs. This model is much less complex than the MS model, but the corresponding bodies have similar length, mass and inertia. We compared the magnitude of torques associated to the human a n d t h e r o b o t movements. The effort distribution on the lower body joints and the torque trajectories of the upper body joints showed interesting similarities. For example, as shown in figure 3, the joint torques exerted at the torso of the robot present the same shape and variation as the human torso torques. However, an important difference is observed and concerns the torque variation along the movement. Indeed, in the human motion, the torques v a r y rather smoothly whereas the robot torques trajectories present sudden peaks at the beginning of each phase of the movement. This is mainly due to the choice of the task dynamics in the sequencing of tasks. Indeed, to guarantee the convergence of each phase of the movement we imposed an exponential decay of the task-functions as it is usually done in robotics. This leads to a succession of exponential decays of torques with a high initial value a t t h e b e g i n n i n g o f e a c h p h a s e . To the contrary, despite the important change in the successive target positions, the human torques appear rather smooth, showing the importance of the anticipation process in human motor control.

Conclusion

This study shows that our advanced motion generation software allows to generate complex movements resembling those of humans. The analysis of the reconstructed human movement and the simulated robot movement shows that the torque variation of the human movement is smoother. This important point suggests that the control design in the task space should be improved in order to take into account the dynamics of whole movement, instead of regulating each task individually.

  The human movement recordings and t h e r o b o t motion generation were performed for t h e s a m e reference movement. The subject was asked to execute a sequence of three reaching tasks with his right hand while b a l a n c i n g o n i t s r i g h t f o o t . Additionally, he was holding a cup in his left hand and must maintain it around a fixed position. The first target was located on the ground in front of the subject. The second target was located at midheight on his right side, and the third one was behind him at shoulder level. Implicitly, d u r i n g each reaching task the subject was required to look at the target. Four snapshots of the successive motion phases are shown in figure1.

Figure 1 :

 1 Figure 1: Snapshots of the motion of the human and the robot model

Figure 2 :

 2 Figure 2: Stack of tasks chronology with task priority decreasing from top to bottom

Figure 3 :

 3 Figure 3: The torques variation of the robot and human torso

Acknowledgments

The authors would like to thank Ko Ayusawa, Yosuke Ikegami and Manish Sreenivasa for making the experiments at the YNL.