Yves Mocquard
email: yves.mocquard@irisa.fr

Bruno Sericola
email: bruno.sericola@inria.fr

Samantha Robert
email: samantha.robert@hotmail.fr

Emmanuelle Anceaume
email: emmanuelle.anceaume@irisa.fr

Analysis of the Propagation Time of a Rumour in Large-scale Distributed Systems

Keywords: rumor spreading, pairwise interactions, Markov chain, analytical performance evaluation. I

The context of this work is the well studied dissemination of information in large scale distributed networks through pairwise interactions. This problem, originally called rumor mongering, and then rumor spreading has mainly been investigated in the synchronous model. This model relies on the assumption that all the nodes of the network act in synchrony, that is, at each round of the protocol, each node is allowed to contact a random neighbor. In this paper, we drop this assumption under the argument that it is not realistic in large scale systems. We thus consider the asynchronous variant, where at time unit, a single node interacts with a randomly chosen neighbor. We perform a thorough study of Tn the total number of interactions needed for all the n nodes of the network to discover the rumor. While most of the existing results involve huge constants that do not allow for comparing different protocols, we prove that in a complete graph of size n ≥ 2, the probability that Tn > k for all k ≥ 1 is less than 1 + 2k(n-2) 2 n 1 -2 n (k-1) . We also study the behavior of the complementary distribution of Tn at point cE(Tn) when n tends to infinity for c = 1. We end our analysis by conjecturing that when n tends to infinity, Tn > E(Tn) with probability close to 0.4484.

important questions of these protocols is the spreading time, that is time it needs for the rumor to be known by all the vertices of the network.

Several models have been considered to answer this question. The most studied one is the synchronous push/pull model, also called the synchronous random phone call model. This model assumes that all the vertices of the network act in synchrony, which allows the algorithms designed in this model to divide time in synchronized rounds. During each synchronized round, each vertex i of the network selects at random one of its neighbor j and either sends to j the rumor if i knows it (push operation) or gets the rumor from j if j knows the rumor (pull operation). In the synchronous model, the spread time of a rumor is defined as the number of synchronous rounds necessary for all the nodes to know the rumor. In one of the first papers dealing with the push operation only, Frieze and Grimmet [START_REF] Frieze | The shortest-path problem for graphs with random arc-lengths[END_REF] proved that if the underlying graph is complete, then asymptotically almost surely the number of rounds is log 2 (n) + log(n) + o(log n) where n is the number of nodes of the graph. Further results have been established (see for example [START_REF] Karp | Randomized rumor spreading[END_REF], [START_REF] Bertier | On gossip and populations[END_REF] and the references herein), the most recent ones resulting from the observation that the rumor spreading time is closely related to the conductance of the graph of the network [START_REF] Giakkoupis | Tight bounds for rumor spreading in graphs of a given conductance[END_REF], [START_REF] Giakkoupis | Tight bounds for rumor spreading with vertex expansion[END_REF]. Investigations have also been done in different topologies of the network [START_REF] Chierichetti | Rumor spreading in social networks[END_REF], [START_REF] Daum | Rumor spreading with bounded indegree[END_REF], [START_REF] Fountoulakis | Rumor spreading on random regular graphs and expanders[END_REF], [START_REF] Panagiotou | Randomized rumor spreading: the effect of the network topology[END_REF], in the presence of link or vertices failures (see [START_REF] Feige | Randomized broadcast in networks[END_REF]), and dynamic graphs [START_REF] Clementi | Rumor spreading in random evolving graphs[END_REF].

All the above studies assume that all vertices of the network act synchronously. In distributed networks, and in particular in large scale distributed systems, such a strict synchronization is unrealistic. Several authors have recently dropped this assumption by considering an asynchronous model. Boyd et al [START_REF] Stephen | Randomized gossip algorithms[END_REF] consider that each node has a clock that goes off at the time of a rate 1 Poisson process. Each time the ring of a node goes off, the push or pull operations are triggered according to the knowledge of the node. Acan et al. [START_REF] Huseyin Acan | On the push&pull protocol for rumour spreading[END_REF] go a step further by studying rumor spreading time for any graph topology. They show that both the average and guaranteed spreading time are Ω(n), where n is the number of nodes in the network. Further investigations have been made for different network topologies [START_REF] Panagiotou | Asynchronous rumor spreading on random graphs[END_REF], [START_REF] Fountoulakis | Ultra-fast rumor spreading in social networks[END_REF].

978-1-5090-3216-7/16/$31.00 c 2016 IEEE a) Our contributions : In this paper we consider the population protocol model, which turns out to resemble to the discrete-time version of the asynchronous spreading model. This model provides minimalist assumptions on the computational power of the nodes: nodes are finite-state automata, identically programmed, they have no identity, they do not know how numerous they are, and they progress in their computation through random pairwise interactions. Their objective is to ultimately converge to a state from which the sought property can be derived from any node [START_REF] Aspnes | An introduction to population protocols[END_REF]. In this model, the spreading time is defined as the number of interactions needed for all the nodes of the network to learn the rumor. Angluin et al [START_REF] Angluin | Fast computation by population protocols with a leader[END_REF] analyze the spreading time of a rumor by only considering the push operation (which they call the oneway epidemic operation), and show that with high probability, a rumor injected at some node requires O(n log n) interactions to be spread to all the nodes of the network.

In the present paper we go a step further by considering a more general problem namely, that is all the nodes of the network initially receive an input value, and the objective for each node is to learn the maximal value initially received by any node. Note that the rumor spreading problem is a particular instance of this problem when two input values 1 and 0 are considered respectively representing the knowledge and the absence of knowledge of the rumor. We present a thorough analysis of the number of interactions needed for all the nodes to converge to the correct response. Specifically, we study the expectation, variance and an exact formulation of the distribution of the number of interactions needed to propagate a rumor. This formulation being hardly usable in practice once n becomes too large, a tight bound is derived. This bound is all the more interesting as usual probabilistic inequalities fail to provide relevant results in this case. Finally, we study the asymptotic behavior of the spreading time when the size of the network tends to infinity.

b) Road map: The remainder of this paper is organized as follows. Section II presents the population protocol model. Section III specifies the problem addressed in this work. Analysis of the spreading time is proposed in Section IV, while we study in Section V its asymptotic behavior. We have simulated our protocol to illustrate our theoretical analysis. Finally, Section VI concludes.

II. POPULATION PROTOCOLS MODEL

In this section, we present the population protocol model, introduced by Angluin et al. [START_REF] Angluin | Computation in networks of passively mobile finite-state sensors[END_REF]. This model describes the behavior of a collection of nodes that interact pairwise. The following definition is from Angluin et al [START_REF] Angluin | The computational power of population protocols[END_REF]. A population protocol is characterized by a 6-tuple (Q, Σ, Y, ι, ω, f ), over a complete interaction graph linking the set of n nodes, where Q is a finite set of states, Σ is a finite set of input symbols, Y is a finite set of output symbols, ι : Σ → Q is the input function that determines the initial state of a node, ω : Q → Y is the output function that determines the output symbol of a node, and f : Q × Q → Q × Q is the transition function that describes how two nodes interact and update theirs states. Initially all the nodes start with a initial symbol from Σ, and upon interactions with nodes update their state according to the transition function f . Interactions between nodes are orchestrated by a random scheduler: at each discrete time, any two nodes are randomly chosen to interact with a given distribution. Note that its is assumed that the random scheduler is fair, which means that the interactions distribution is such that any possible interaction cannot be avoided forever. The notion of time in population protocols refers to as the successive steps at which interactions occur, while the parallel time refers to as the successive number of steps each node executes [START_REF] Aspnes | An introduction to population protocols[END_REF]. Nodes do not maintain nor use identifiers (nodes are anonymous and cannot determine whether any two interactions have occurred with the same agents or not). However, for ease of presentation the nodes are numbered 1, 2, . . . , n. We denote by C (i) t the state of node i at time t. The stochastic process C = {C t , t ≥ 0}, where

C t = (C (1) t , . . . , C (n) t )
, represents the evolution of the population protocol. The state space of C is thus Q n and a state of this process is also called a protocol configuration.

III. SPREADING THE MAXIMUM

We consider in this section the following problem. Each site has initially an integer value. At each discrete instant of time, two distinct nodes are successively chosen and they change their value with the maximum value of each node. More precisely, for all nodes a and b, with a = b, we consider the function f given by f (a, b) = (max{a, b}, max{a, b}) .

We want to evaluate the time needed so that all the nodes get the same value.

Let C = {C t , t ≥ 0} be a discrete-time stochastic process with state space S = N n . For every t ≥ 0, the state at time t of the process is denoted by

C t = (C (1) t , . . . , C (n) t ), where C (i)
t is the integer value of node i at time t. At each instant t, two distinct indexes i and j are successively chosen among the set of nodes 1, . . . , n randomly. We denote by X t the random variable representing this choice and we suppose that this choice is uniform, i.e we suppose that

P{X t = (i, j)} = 1 n(n -1) 1 {i =j} .
Once the couple (i, j) is chosen at time t, the process reaches state C t+1 , at time t + 1, given by

C (i) t+1 = C (j) t+1 = max{C (i) t , C (j) t } and C (m) t+1 = C (m) t for i = j.
We denote by M the maximum initial value among all the nodes, i.e. M = max{C

(1) 0 , . . . , C (n) 0 }. It is easily checked that for all t ≥ 0, we have M = max{C (1) t , . . . , C (n) t }.
We consider the random variable T n defined by

T n = inf{t ≥ 0 | C (i) t = M, for every 1, . . . , n}.
The random variable T n represents the number of interactions needed for all the nodes in the network to know the maximal value M . We introduce the discrete-time stochastic process Y = {Y t , t ≥ 0} with state space {1, . . . , n} defined, for all t ≥ 0, by

Y t = i | C (i) t = M .
The random variable Y t represents the number of nodes knowing the maximum value M at time t. The stochastic process Y is then a homogeneous Markov chain with transition probability matrix A. The non zero transition probabilities are given, for i, j = 1, . . . , n, by

         A i,i = 1 - 2i(n -i) n(n -1) , A i,i+1 = 2i(n -i) n(n -1) , for i = n.
Indeed, when Y t = i, in order to get Y t+1 = i + 1, either the first node must be chosen among the ones with the maximum value (probability i/n) and the second agent must be chosen among the ones with the non maximum value (probability (ni)/(n -1)) or the first node must be chosen among the ones with the non maximum value (probability (n -i)/n) and the second node must be chosen among the ones with the non maximum value (probability i/(n -1)).

The states 1, . . . , n -1 of Y are transient and state n is absorbing. The random variable T n can then be written as

T n = inf{t ≥ 0 | Y t = n}.
It is well-known, see for instance [START_REF] Sericola | Markov Chains. Theory, Algorithms and Applications[END_REF], that the distribution of T n is given, for every k ≥ 0, by

P{T n > k} = αQ k 1, ( 1
)
where α is the row vector containing the initial probabilities of states 1, . . . , n -1, that is α i = P{Y 0 = i}, Q is the submatrix obtained from A by deleting the row and the column corresponding to absorbing state n and 1 is the column vector of dimension n -1 with all its entries equal to 1.

For i = 0, . . . , n, we introduce the notation

p i = 2i(n -i) n(n -1)
and we denote by H k the harmonic series defined by H 0 = 0 and

H k = k =1 1/ , for k ≥ 1.
Note that, for every i = 0, . . . , n, we have p i = p n-i .

If we denote by S i , for i = 1, . . . , n -1, the total time spent by the Markov chain Y in state i, then conditionally on the event Y 0 = i, S has a geometric distribution with parameter p , for = i, . . . , n -1 and in this case, we have

T n = S i + • • • + S n-1 . It follows that P{T n > k | Y 0 = i} = P{S i + • • • + S n-1 > k}, which means that P{T n > k | Y 0 = i} is decreasing with i and in particular that P{T n > k | Y 0 = i} ≤ P{T n > k | Y 0 = 1}.
(2)

IV. ANALYSIS OF THE SPREADING TIME

In the following we study the expectation and the variance of T n , the number of interactions needed for all the nodes in the network to know the maximal value M . We then provide an explicit expression of the distribution of T n , and then a bound and an equivalent for the explicit distribution of T n .

A. Expectation and variance of T n

The mean time E(T n ) needed so that all the nodes get the same value is then given by

E(T n ) = α(I -Q) -1 1, ( 3
)
where I is the identity matrix. This expectation can also be written as

E(T n ) = n-1 i=1 α i E(T n | Y 0 = i).
This conditional expectations are given by the following theorem.

Theorem 1: For every n ≥ 1 and i = 1, . . . , n, we have

E(T n | Y 0 = i) = (n -1)(H n-1 + H n-i -H i-1 )
2 .

Proof: If Y 0 = n, which means that all the nodes start with same values, then we have

T n = 0 and so E(T n | Y 0 = n) = 0. For i = 1, . . . , n -1 we have E(T n | Y 0 = i) = n-1 =i E(S ) = n-1 =i 1 p = n(n -1) 2 n-1 =i 1 (n -) = n -1 2 n-1 =i 1 + 1 n - = (n -1)(H n-1 + H n-i -H i-1 ) 2 ,
which completes the proof.

In particular, when the maximum value is initially unique, i.e. when Y 0 = 1 with probability 1, we have α 1 = 1 and thus

E(T n ) = E(T n | Y 0 = 1) = (n -1)H n-1 ∼ n-→∞ n ln(n).
More generally, from Relation (2), we have

E(T n ) ≤ E(T n | Y 0 = 1) = (n -1)H n-1 ∼ n-→∞ n ln(n).
The variance of T n is obtained similarly.

Theorem 2: For every n ≥ 1 and i = 1, . . . , n, we have

Var(T n | Y 0 = i) = (n -1) 2 4 n-1 =i 1 2 + n-i =1 1 2 - E(T n | Y 0 = i) n .
Proof: If Y 0 = n, which means that all the nodes start with the same values, then we have T n = 0 and thus Var(T n | Y 0 = n) = 0. For i = 1, . . . , n -1 we have, using the independence of the S ,

Var(T n | Y 0 = i) = n-1 =i Var(S ) = n-1 =i 1 -p p 2 = n-1 =i 1 p 2 - n-1 =i 1 p = n 2 (n -1) 2 4 n-1 =i 1 2 (n -) 2 - n(n -1) 2 n-1 =i 1 (n -) = (n -1) 2 4 n-1 =i 1 + 1 n - 2 - n(n -1) 2 n-1 =i 1 (n -) = (n -1) 2 4 n-1 =i 1 2 + 1 (n -) 2 - n -1 2 n-1 =i 1 (n -) = (n -1) 2 4 n-1 =i 1 2 + 1 (n -) 2 - E(T n | Y 0 = i) n = (n -1) 2 4 n-1 =i 1 2 + n-i =1 1 2 - E(T n | Y 0 = i) n ,
which completes the proof.

In particular, when the maximum value is initially unique, i.e. when Y 0 = 1 with probability 1, we have α 1 = 1 and thus

Var(T n ) = Var(T n | Y 0 = 1) = (n -1) 2 2 n-1 =1 1 2 - n -1 n H n-1 ∼ n-→∞ π 2 n 2 12 .
More generally, from Theorem 2, we have

Var(T n | Y 0 = i) ≤ (n -1) 2 4 n-1 =i 1 2 + n-i =1 1 2 ≤ (n -1) 2 2 n-1 =1 1 2 ≤ π 2 n 2 12 .

It follows that

Var(T n ) = n-1 i=1 α i Var(T n | Y 0 = i) ≤ π 2 n 2 12 .

B. Explicit expression of the distribution of T n

The distribution of T n , for n ≥ 2, which is given by Relation (1) can be easily computed as follows. Let

V (k) = (V 1 (k), . . . , V n-1 (k)) be the column vector defined by V i (k) = P{T n > k | Y 0 = i}. According to Rela- tion (1), we have V (k) = Q k 1. Since V (0) = 1, writing V (k) = QV (k -1) for k ≥ 1, we get for i = 1, . . . , n -2, V i (k) = (1 -p i ) V i (k -1) + p i V i+1 (k -1), V n-1 (k) = (1 -p n-1 ) V n-1 (k -1). (4)
Recall that we have p i = 2i(n -i)/(n(n -1)). This recursion can be easily computed since we have, for k ≥ 0,

V n-1 (k) = (1 -p n-1 ) k = 1 - 2 n k . (5)
In the next theorem, we derive from recursion (4) an explicit expression of the distribution of T n .

Theorem 3: For every n ≥ 1, k ≥ 0 and i = 1, . . . , n -1, we have

P{T n > k |Y 0 = n -i} = n/2 j=1 (c i,j (1 -p j ) + kd i,j ) (1 -p j ) k-1 ,
where the coefficients c i,j and d i,j , which do not depend on k, are given, for j = 1, . . . , n -1, by c 1,j = 1 {j=1} and d 1,j = 0 and for i ∈ {2, . . . , n -1} by

c i,j = p i c i-1,j p i -p j - p i d i-1,j (p i -p j ) 2 for i = j, n -j, d i,j = p i d i-1,j p i -p j for i = j, n -j, c i,i = 1 - n/2 j=1,j =i c i,j for i ≤ n/2, c i,n-i = 1 - n/2 j=1,j =n-i c i,j for i > n/2, d i,i = p i c i-1,i for i ≤ n/2, d i,n-i = p i c i-1,n-i for i > n/2.
Proof: See appendix

C. Bounds of the distribution of T n

The exact expression of the distribution of T n presented earlier is hardly usable in practice, and computation using formula (4) may take a long time for large values of n. To overcome this problem, we propose in this section a bound and an equivalent for the quantity P{T n > k | Y 0 = i} derived from the recursive formula (4).

Theorem 4: For all n ≥ 2 and k ≥ 1 we have

P{T n > k | Y 0 = 1} ≤ 1 + 2k(n -2) 2 n 1 - 2 n k-1 , P{T n > k | Y 0 = 1} ∼ k-→∞ 1 + 2k(n -2) 2 n 1 - 2 n k-1
and for i = 2, . . . , n -1 and k ≥ 0,

P{T n > k | Y 0 = i} ≤ (n -i)(n -2) i -1 1 - 2 n k , P{T n > k | Y 0 = i} ∼ k-→∞ (n -i)(n -2) i -1 1 - 2 n k .
Moreover, we have

P{T n > k} ≤ P{T n > k | Y 0 = 1}.
Proof: The result is trivial for n = 2 since in this case we have T 2 = 1. We thus suppose that n ≥ 3. Note that by definition of p i we have p i = p n-i . Consider the sequence b i defined for i = 1, . . . , n -2, by

b 1 = 1 and b i = p i b i-1 p i -p 1 , for i = 2, . . . , n -2.
Observing that

b i = i(n -i)b i-1 (i -1)(n -i -1)
, it is easily checked by recurrence that for i = 1, . . . , n -2, we have

b i = i(n -2) n -i -1 .
We show now by recurrence that for all i = 1, . . . , n -2, we have

V n-i (k) ≤ b i (1 -p 1 ) k , for all k ≥ 0 and V n-i (k) ∼ k-→∞ b i (1 -p 1 ) k .

Both results are true for

i = 1 since V n-1 (k) = (1-p n-1 ) k = (1 -p 1 )
k . Suppose now that these results are true for a fixed integer i with 1 ≤ i ≤ n -3. From Relations (4), we have

V n-i-1 (k) = (1 -p n-i-1 )V n-i-1 (k -1) + p n-i-1 V n-i (k -1) = (1 -p i+1 )V n-i-1 (k -1) + p i+1 V n-i (k -1).
Using the recurrence hypothesis, we obtain, for what concerns the inequality,

V n-i-1 (k) ≤ (1 -p i+1 )V n-i-1 (k -1) + p i+1 b i (1 -p 1 ) k-1 .
Expanding this inequality and using the fact that V n-i-1 (0) = 1, this leads to

V n-i-1 (k) ≤ (1 -p i+1 ) k + p i+1 b i k-1 j=0 (1 -p i+1 ) j (1 -p 1 ) k-1-j = (1 -p i+1 ) k + p i+1 b i (1 -p 1 ) k -(1 -p i+1 ) k p i+1 -p 1 = (1 -p i+1 ) k + b i+1 (1 -p 1 ) k -(1 -p i+1 ) k (1 -b i+1 ) (1 -p i+1 ) k + b i+1 (1 -p 1 ) k . Since b i+1 ≥ 1, we get V n-i-1 (k) ≤ b i+1 (1 -p 1 ) k
In the same way, using a similar calculus, we obtain

V n-i-1 (k) ∼ k-→∞ (1 -b i+1 ) (1 -p i+1 ) k + b i+1 (1 -p 1 ) k .
Since p i+1 > p 1 , we also get

V n-i-1 (k) ∼ k-→∞ b i+1 (1 -p 1 ) k .
We thus have shown that for all i = 1, . . . , n -2, we have

V n-i (k) ≤ b i (1 -p 1 ) k , for all k ≥ 0 and V n-i (k) ∼ k-→∞ b i (1 -p 1 ) k .
In particular, for i = n -2 we obtain

V 2 (k) ≤ b n-2 (1 -p 1 ) k , for all k ≥ 0 and V 2 (k) ∼ k-→∞ b n-2 (1 -p 1 ) k .
Consider now the term V 1 (k). From Relations (4) and using the previous inequality, we have

V 1 (k) = (1 -p 1 )V 1 (k -1) + p 1 V 2 (k -1) ≤ (1 -p 1 )V 1 (k -1) + p 1 b n-2 (1 -p 1 ) k-1 .
Expanding this inequality and using the fact that V 1 (0) = 1, this leads to

V 1 (k) ≤ (1 -p 1 ) k + p 1 b n-2 k-1 j=0 (1 -p 1 ) j (1 -p 1 ) k-1-j = (1 -p 1 ) k + p 1 b n-2 k(1 -p 1 ) k-1 = (1 -p 1 + kp 1 b n-2 ) (1 -p 1 ) k-1 ≤ (1 + kp 1 b n-2 ) (1 -p 1 ) k-1 , which gives V 1 (k) ≤ 1 + 2k(n -2) 2 n 1 - 2 n k-1 .
In the same way, using a similar calculus, we obtain

V 1 (k) ∼ k-→∞ 1 + 2k(n -2) 2 n 1 - 2 n k-1 . Finally, since P{T n > k | Y 0 = i} is decreasing with i, we have P{T n > k} = n-1 i=1 P{T n > k | Y 0 = i}P{Y 0 = i} ≤ P{T n > k | Y 0 = 1},
which completes the proof.

The bound established in Theorem 4 is all the more interesting as usual probabilistic inequalities fail to provide relevant results in this particular case. For example, Markov inequality leads for all real number c ≥ 1 to

P{T n ≥ cE(T n )} ≤ 1 c ,
and Bienaymé-Tchebychev inequality leads for all real number x > 0 to

P{|T n -E(T n )| ≥ x} ≤ π 2 n 2 12x 2 .
The author of [START_REF] Janson | Tail bounds for sums of geometric and exponential variables[END_REF] provides a bound, based on Chernoff inequality, for the tail probabilities of the sum of independent, but not necessarily identically distributed, geometric random variables. In the particular case of our protocol computing the maximum, this leads to the following result.

Theorem 5: For all n ≥ 3 and for all real number c ≥ 1, we have

P(T n > cE(T n )) ≤ 1 c 1 - 2 n (c -1 -ln c)(n -1)H n-1 .
The right-hand side term is equal to 1 when c = 1.

Proof: We have already shown that

P(T n > cE(T n )) ≤ P(T n > cE(T n ) | Y 0 = 1}.
The upper bound is then an application of Theorem 2.3 of [START_REF] Janson | Tail bounds for sums of geometric and exponential variables[END_REF], and it is clearly equal to 1 when c = 1.

Applying Theorem 4 at point k = cE(T n ) , we obtain

P(T n > cE(T n )) ≤ 1 + 2 cE(T n ) (n -2) 2 n × 1 - 2 n cE(Tn) -1 ≤ 1 + 2cE(T n )(n -2) 2 n × 1 - 2 n cE(Tn)-2
.

From now on we denote this bound by f (c, n) and in the same way, we denote by g(c, n) the bound of P(T n > cE(T n ))

derived from Theorem 5. We then have, for n ≥ 3 and c ≥ 1,

f (c, n) = 1 + 2c(n -1)H n-1 (n -2) 2 n × 1 - 2 n c(n-1)Hn-1-2 g(c, n) = 1 c 1 - 2 n (c-1-ln(c))(n-1)Hn-1
.

We also introduce the notation

e(c, n) = P(T n > cE(T n )).
Theorem 6: For every n ≥ 3, there exists a unique c * ≥ 1 such that f (c * , n) = g(c * , n) and we have

f (c, n) > g(c, n) for all 1 ≤ c < c * f (c, n) < g(c, n) for all c > c * . (6) Furthermore, lim c-→∞ f (c, n) g(c, n) = 0.
Proof: See appendix.

The graphs on Figures 1, 2 and 3 illustrate the behavior of the bounds f (c, n) and g(c, n), depending on c and for different values of n, compared to the real distribution of T n at point cE(T n ), i.e. to e(c, n) = P{T n > E(T n )}. The bound f (c, n) that we provided in Theorem 4 clearly shows better accuracy than the Chernoff bound g(c, n) provided in [START_REF] Janson | Tail bounds for sums of geometric and exponential variables[END_REF] above the threshold c * introduced in Theorem 6. Furthermore, this threshold seems to decrease to 1 as n tends to infinity, as can be seen on Figure 4.

V. ASYMPTOTIC ANALYSIS OF THE DISTRIBUTION OF T n

We analyze in this section the behavior of the complementary distribution of T n at point cE(T n ) when n tends to infinity, depending on the value of c.

We prove in the following corollary that the bounds f (c, n) and g(c, n), obtained from Theorem 4 and Theorem 5 respectively with k = cE(T n ), both tend to 0 when n goes to infinity. Proof: For all x ∈ [0, 1), we have ln(1 -x) ≤ -x. Applying this property to the bound f (c, n) leads to

f (c, n) ≤ 1 + 2c(n -1)H n-1 (n -2) 2 n ×e -2(c(n-1)Hn-1-2)/n ≤ 1 + 2c(n -2) 2 H n-1 e -2(c(n-1)Hn-1-2)/n . Since ln(n) ≤ H n-1 ≤ 1 + ln(n -1), we get f (c, n) ≤ 1 + 2c(n -2) 2 (1 + ln(n -1) ×e -2(c(n-1) ln(n)-2)/n = 1 + 2c(n -2) 2 (1 + ln(n -1)) e -2c ln(n) ×e 2(c ln(n)+2)/n .
For x ≥ 0, the function u(x) = e 2(c ln(x)+2)/x satisfies u(x) ≤ exp 2c/e (c-2)/c , so we obtain

f (c, n) ≤ 1 + 2c(n -2) 2 (1 + ln(n -1)) n 2c exp 2c/e (c-2)/c .
The fact that c > 1 implies that this last term tends to 0 when n -→ ∞. Concerning the bound g(c, n), we have

g(c, n) = 1 c 1 - 2 n (c-1-ln(c))(n-1)Hn-1 = 1 c e (c-1-ln(c))(n-1)Hn-1 ln(1-2/n) ≤ 1 c
e -2(c-1-ln(c))(n-1)Hn-1/n , which tends to 0 when n tends to infinity, since c -1 -ln(c) is positive for c > 1.

Theorem 8: For all real c ≥ 0, we have

lim n→+∞ P{T n > cE(T n )} = 0 if c > 1 1 if c < 1.
Proof: From Corollary 7, both bounds f (c, n) and g(c, n)

of P{T n > cE(T n )} tend to 0 when n tends to infinity, so using either f (c, n) or g(c, n) we deduce that lim n-→∞ P{T n > cE(T n )} = 0 for all c > 1.

In the case where c < 1, Theorem 3.1 of [START_REF] Janson | Tail bounds for sums of geometric and exponential variables[END_REF] leads to

P{T n > cE(T n )} ≥ 1 -e -2(n-1)Hn-1(c-1-ln(c))/n ≥ 1 -e -2(n-1) ln(n)(c-1-ln(c))/n .
Since c -1 -ln(c) > 0 for all c ∈ [0, 1), the right-hand side term of this inequality tends to 1 when n -→ ∞. Thus, In Figure 5, we show the probability P{T n > E(T n )} for different values of n. The oscillations of this probability with n are due to the fact T n is a discrete random variable and E(T n ) is not an integer. That is why we propose in this figure a smoothing of this probability using the sequence

lim n-→∞ P{T n > cE(T n )} = 1 when c < 1.
P{T n > E(T n )} (1-a n )P{T n > E(T n )} + a n P{T n > E(T n )+1} P{T n > E(T n )+1}
s n = (1 -a n )P{T n > E(T n )} + a n P{T n > E(T n ) + 1}
, where a n is the fractional part of E(T n ), that is a n = E(T n )-E(T n ) . Since a n ∈ [0, 1], we have

P{T n > E(T n ) + 1} ≤ s n ≤ P{T n > E(T n )},
that is why we also show in this figure the probability P{T n > E(T n ) + 1}. We also checked that the sequence (s n ) is increasing until n = 20000. This figure suggests the following result proposed as a conjecture.

Conjecture : lim n-→∞ P{T n > E(T n )} exists and ≈ 0.4484.

Figure 6 shows the results obtained by simulation concerning the proportion of nodes informed by rumor as a function of the parallel time. Recall that the parallel time refers to as the successive number of steps each node executes [START_REF] Aspnes | An introduction to population protocols[END_REF]. Initially, a single node is informed of the rumor. This figure illustrates our analysis. For instance, with probability almost 1 one thousand nodes (resp. one million nodes) learn the rumor after no more than 7 (resp 11) interactions for each of them. The complexity in space (number of memory bits) is in O(1).

VI. CONCLUSION

In this paper we have provided a thorough analysis of the rumor spreading time in the population protocol model. Providing such a precise analysis is a step towards the design of more complex functionality achieved by combining simple population protocols [START_REF] Michail | Terminating population protocols via some minimal global knowledge assumptions[END_REF], [START_REF] Angluin | Fast computation by population protocols with a leader[END_REF]. Indeed, an important feature of population protocols is that they do not halt. Nodes can never know whether their computation is completed and thus nodes forever interact with their neighbors while their outputs stabilize to the desired value (e.g. the maximal value of any node of the network). By precisely characterizing, for each protocol of interest, with any high probability, the number of interactions each node must execute to converge to the desired value, each node can on its own, decide the time from which the current protocol has stabilized and start the parallel of sequential executions of the next ones.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Bounds f (c, n) and g(c, n) beside the real value of P(Tn > cE(Tn)) = e(c, n) for n = 100, as functions of c. In this case, we have c * = 1.14641.

Fig. 3 .

 3 Fig. 3. Bounds f (c, n) and g(c, n) beside the real value of P(Tn > cE(Tn)) = e(c, n) for n = 5000, as functions of c. In this case, we have c * = 1.11385.

7 cFig. 4 .

 74 Fig. 4. Approximate values of c * for different network sizes n.

Corollary 7 :

 7 For all real number c > 1, we havelim n-→∞ f (c, n) = 0 and lim n-→∞ g(c, n) = 0.

Fig. 5 .Fig. 6 .

 56 Fig. 5. P{Tn > E(Tn)} as a function of n and its smoothing obtained with an = E(Tn) -E(Tn) .

One of the This work was partially funded by the French ANR project SocioPlug (ANR-13-INFR-0003), and by the DeSceNt project granted by the Labex CominLabs excellence laboratory (ANR-10-LABX-07-01).

APPENDIX

We give in this appendix the proofs of Theorem 3 and Theorem 6. In order to prove Theorem 3, we first need the following Lemma.

Lemma 9: Let N ≥ 1, a ∈ (0, 1), b 1 , . . . , b N ∈ (0, 1), c 1 , . . . , c N ∈ R and d 1 , . . . , d N ∈ R, with the condition, for every j = 1, . . . , N, d j = 0 if b j = a.

Then the sequence (u k ) k≥0 defined by u 0 = 1 and u k+1 = au k + N j=1 (c j b j + kd j )b k-1 j , k ≥ 0 [START_REF] Bertier | On gossip and populations[END_REF] satisfies

where

Proof: We prove this lemma by recurrence. For k = 0, Relation (8) gives u 0 = 1. We introduce the notation α = 1 -N j=1 θ j 1 {bj =a} and f j (k) = θ j b j + kγ j . Relation [START_REF] Censor-Hillel | Global computation in a poorly connected world: Fast rumor spreading with no dependence on conductance[END_REF] can then be rewritten as

Suppose that this last relation is true for a fixed k ≥ 0. From Relation [START_REF] Bertier | On gossip and populations[END_REF], we obtain

Writing c j b j = c j b j 1 {bj =a} +c j b j 1 {bj =a} and d j = d j 1 {bj =a} since d j = 0 when b j = a, we obtain

The first term of this last summation can be simplified as follows. By definition of f j (k) and observing that aγ j + d j = γ j b j , we have

= aθ j b j + c j b j + kγ j b j = (aθ j + c j + kγ j )b j = (aθ j + c j -γ j + (k + 1)γ j )b j .

Since c j -γ j = (b j -a)θ j , this last expression leads to

Putting this result into [START_REF] Chierichetti | Rumor spreading in social networks[END_REF] gives

We are now ready to prove Theorem 3.

where the coefficients c i,j and d i,j , which do not depend on k, are given, for j = 1, . . . , n -1, by c 1,j = 1 {j=1} and d 1,j = 0 and for i ∈ {2, . . . , n -1} by

Proof of Theorem 3: The proof is made by recurrence on integer i. In fact, we prove that for every i = 1, . . . , n -1, we have

Relations ( 11) and ( 12) are true for i = 1 since by definition we have c 1,j = 1 {j=1} and d 1,j = 0. It follows that Relation [START_REF] Clementi | Rumor spreading in random evolving graphs[END_REF] gives, for i = 1, V n-1 (k) = (1 -p 1 ) k , which is in accordance with Relation (4).

Suppose now that Relations [START_REF] Clementi | Rumor spreading in random evolving graphs[END_REF], [START_REF] Daum | Rumor spreading with bounded indegree[END_REF] and ( 12) are true for a fixed integer i, i ≤ n -2. Using Relation (4) at point k + 1 and the fact that p n-i-1 = p i+1 , we obtain

Integer i being fixed, we apply Lemma (9) in which we set u k = V n-i-1 (k), a = 1 -p i+1 , N = n/2 and for j = 1, . . . , N , b j = 1 -p j , c j = p i+1 c i,j and d j = p i+1 d i,j . Note that condition of Lemma ( 9) is satisfied. Indeed, a = b j is equivalent to either i + 1 = j or i + 1 = n -j. Since i + 1 = j is equivalent to i + j + 1 = 2j and since j ≤ n/2, both conditions implies that i + j < n which means from Relation (11) that d j = p i+1 d i,j = 0. With this notation, the parameters θ j and γ j writes

Applying Lemma ( 9) from ( 13), we obtain

we obtain, using again the fact that p i+1 = p n-i-1 ,

Now, we must prove the recurrence for the two additional properties [START_REF] Daum | Rumor spreading with bounded indegree[END_REF] and [START_REF] Demers | Epidemic algorithms for replicated datbase maintenance[END_REF]. Note that i + 1 < j implies i + 1 = j and i + 1 = n -j. So if i + 1 < j then we have i < j and i < n/2 < n -j which means that c i,j = d i,j = 0, which in turn implies that θ j = 0 and thus c i+1,j = 0

In the same way, note that i+1 < n-j implies i+1 = n-j. So if i + 1 < n -j then we have i < n -j which means that d i,j = 0, which in turn implies that γ j = 0 and thus d i+1,j = 0. This completes the proof.

Theorem 6 For every n ≥ 3, there exists a unique c * ≥ 1 such that f (c * , n) = g(c * , n) and we have

Proof of Theorem 6: In order to simplify the writing, we introduce the following notations.

Firstly, since

Taking the derivative with respect to c, gives

The term C n c An is strictly positive for all c ≥ 1, so we have

It is easy to check that -(A n + 1)/(B n (A n + 2)) < 0 for all n ≥ 3, which means that ∂(f /g)/∂c < 0 for all c ≥ 1. This implies that the function c -→ f (c, n)/g(c, n) is strictly decreasing on [1, +∞).

Observing that A n < -2 for every n ≥ 3, we obtain

which proves the second part of the theorem. Secondly, since

it is easy to check that f (1, n)/g(1, n) ≥ 1, for all n ≥ 3.

Let us recapitulate and conclude. For all n ≥ 3, we have f (1, n)/g(1, n) ≥ 1 and c -→ f (c, n)/g(c, n) is continuous and strictly decreasing to 0. It follows that there exists a unique value c * ≥ 1 such that f (c * , n)/g(c * , n) = 1 and satisfying the conditions [START_REF] Fountoulakis | Rumor spreading on random regular graphs and expanders[END_REF]. This completes the proof.


François Vignale

EuRED : XML-TEI in database structure

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

  EuRED : XML-TEI in database structure The EuRED collects testimonies of European Reading Experiences from the 15th century until today through text mining, public and scholarly crowdsourcing and contribution interface.This multilingual, multimedia, open access and interoperable database, inspired by the 2006 British UK-RED, includes inquiries, geographical and chronological extensions. While addressing new sources (CMCs, audio transcriptions) and interests by the general public, policy makers, educators and scholars, this innovative tool aims to offer access to the evolving journey of reading. Technical context and specifications -Proof of concept -Enrichment of existing data -Interoperability -Linked open data -Open science principles Methods -Definition of what is a reading exprience -Work on a new data model -Choice of CollectiveAccess as CMS -Choice of XML-TEI for data encoding XML-TEI in EuRED -Creation of a new class (model.experienceDescPart) -Creation of 23 new elements mostly describing the reading experience -Description of the reading experience placed into the <teiHeader/> -Use of pointers into <text/> referring to the reading experience What has been done -New data model based on XML-TEI -Creation of specific thesaurii -Conversion of UK-RED data into TEI -Linking with referentials (ISO, Geonames ...) What has to be done -Full disambiguation -Automatic enrichment -Linking with more referentials (ISNI, Dbpedia) -Crowdsourcing process François VIGNALE, 3LAM, Université du Maine, Le Mans, France More about the project : http://eured.univ-lemans.fr http://github.com/EuRED Bibliothèques universitaires Le Mans Laval


Dennis Christie
email: dennischristie@etu.u-bourgogne.fr

Cansen Jiang
email: cansen.jiang@u-bourgogne.fr

Danda Paudel
email: danda-pani.paudel@u-bourgogne.fr

Cédric Demonceaux
email: cedric.demonceaux@u-bourgogne.fr

3D Reconstruction of Dynamic Vehicles using Sparse 3D-Laser-Scanner and 2D Image Fusion

Keywords: 3D Reconstruction, 2D camera, 3D camera, Registration, Point Cloud, ICP

Map building becomes one of the most interesting research topic in computer vision field nowadays. To acquire accurate large 3D scene reconstructions, 3D laser scanners are recently developed and widely used. They produce accurate but sparse 3D point clouds of the environments. However, 3D reconstruction of rigidly moving objects along side with the largescale 3D scene reconstruction is still lack of interest in many researches. To achieve a detailed object-level 3D reconstruction, a single scan of point cloud is insufficient due to their sparsity. For example, traditional Iterative Closest Point (ICP) registration technique or its variances are not accurate and robust enough to registered the point clouds, as they are easily trapped into the local minima. In this paper, we propose an 3-Point RANSAC with ICP refinement algorithm to build 3D reconstruction of rigidly moving objects, such as vehicles, using 2D-3D camera setup. Results show that the proposed algorithm can robustly and accurately registered the sparse 3D point cloud.

I. INTRODUCTION A. Background

Map building becomes one of the most interesting research in computer vision field nowadays, due to its wide future applications such as localization, navigation, and autonomous driving vehicles. Simultaneously Localization and Mapping (SLAM) technique is widely used by researchers to build the 3D map of an environment. Multiple cameras and odometry sensors are involved to serve the system, which allow us to acquire 3D data of the scene to build the 3D map of the environment. A benchmark, provide by Geiger et al. [START_REF] Geiger | Are we ready for autonomous driving? the kitti vision benchmark suite[END_REF], utilized 3D laser scanner as well as 2D cameras on the vehicle, which produces very accurate but sparse 3D measurement of the scene. SLAM approach can be used to build a 3D reconstruction for the static objects such as buildings or streets.

However, in the real world, there are usually not only static objects, but also moving objects in outdoor environments. It is very common that the outdoor environment is a dynamic scene, for instance, people are walking around, vehicles are running, cyclists are crossing the roads, which make the outdoor environment barely a completely static scene. When there are a large number of dynamic objects around, 3D map building by SLAM will not be accurate. 3D reconstruction in a dynamic scene is still an unsolved and challenging task.

A complete framework of static map and dynamic object reconstruction in outdoor environments was proposed by Jiang et al. [START_REF] Jiang | Static-map and dynamic object reconstruction in outdoor scenes using 3d motion segmentation[END_REF], [START_REF] Danda | 2d-3d camera fusion for visual odometry in outdoor environments[END_REF]. They succeeded to detect and segment the dynamic objects from the scene to reconstruct the static scene parts as well as the dynamic objects independently. The proposed static map building algorithm achieved very satisfactory performances on realistic outdoor environment.

This paper is a complementary of [START_REF] Jiang | Static-map and dynamic object reconstruction in outdoor scenes using 3d motion segmentation[END_REF] with a smaller scope. We mainly focus on improving the reconstruction quality of the dynamic objects from a moving camera setup. As the output of the framework, multiple dynamic objects, such as vehicles, are obtained from different view ports due to the relative motions between the camera and the objects. In such kind of scenarios, reconstructing the moving objects from a moving camera setup is very challenging. To tackle this problem, the objective of this work is to build a 3D reconstruction of rigidly moving vehicle, using a 2D-3D calibrated camera setup.

B. Problem Definition

To reconstruct the rigidly moving objects from a moving camera setup, we propose to use 3D laser-scanner, such as Velodyne [START_REF] Moosmann | Velodyne slam[END_REF], due to its stability and accuracy in acquiring the 3D data. Yet, the laser scanner provides only sparse 3D point cloud, especially for the object that is far away from the camera. Moreover, the laser scanner does not give us texture information of the point cloud. We define a high quality 3D reconstruction as: object is reconstructed accurately, densely and textured. The dilemma comes when texture is obtained by using 2D cameras, while the 3D points are acquired from laser scanner. A calibrated 2D-3D camera setup can overcome this dilemma, which means that there exist correspondences between 2D image points with 3D laser points. For instance, giving a mobile robot equipped with calibrated 2D and 3D cameras, one can reconstruct photometric high quality 3D scenes for static objects by registering a sequence of observations [START_REF] Jiang | Static-map and dynamic object reconstruction in outdoor scenes using 3d motion segmentation[END_REF], [START_REF] Danda | 2d-3d camera fusion for visual odometry in outdoor environments[END_REF]. Differently, registering the moving objects are more challenging due to their sparsity, small overlapping area, and non-uniform point density distribution of the moving objects. Two typical failure cases are shown in Fig. 1 and Fig. 2. The registrations using traditional ICP algorithm are failed due to the reasons of poor 3D geometry structure and small partially overlapping, respectively. This is due to traditional ICP algorithm only consider the mean square error (MSE) minimization between points of two point clouds. Minimum MSE, doesn't guarantee it converge into correct registration.

To overcome problem that is stated, we propose a series of registration algorithm, a 3D reconstruction of dynamic objects using 2D-3D camera setup attached to a moving vehicle, as the contribution of this paper.

II. LITERATURE REVIEW Finding feature correspondences is a fundamental problem in computer vision applications, such as stereo vision [START_REF] Marr | A computational theory of human stereo vision[END_REF], struct-from-motion [START_REF] Robert C Bolles | Epipolarplane image analysis: An approach to determining structure from motion[END_REF], motion analysis [START_REF] Michela | Motion analysis by feature tracking[END_REF], and especially point cloud registration [START_REF] Rusinkiewicz | Efficient variants of the icp algorithm[END_REF]. In literature, there are many methods to match the features, such as statistical approach with similarity measurement [START_REF] Candocia | A similarity measure for stereo feature matching[END_REF], feature descriptor based matching [START_REF] Montesinos | Matching color uncalibrated images using differential invariants[END_REF]- [START_REF] David G Lowe | Distinctive image features from scale-invariant keypoints[END_REF].

One of the most robust descriptor based feature matching algorithm is Scale Invariant Feature Transform (SIFT) by Lowe [START_REF] David G Lowe | Distinctive image features from scale-invariant keypoints[END_REF]. SIFT is able to detect keypoints and build the scale invariant and orientation invariant descriptors. Although descriptor based feature matching is robust and efficient in rich-texture regions, features from low-texture area are easily wrongly matched. Thus, feature descriptor matching is generally reliable in matching highly distinctive features.

For feature matching between successive frames, KLT [START_REF] Tomasi | Detection and tracking of point features[END_REF] feature tracking algorithm is applied to establish the dense correspondences [START_REF] Jiang | Static-map and dynamic object reconstruction in outdoor scenes using 3d motion segmentation[END_REF]. The main advantage of OF feature tracking and matching comes from its ability to produce very dense feature matches unlike the sparse feature matches from feature descriptors. In addition, a cross validation can be applied to improve the tracking accuracy and robustness. Since image sequence (video) is involved in this work, optical flow is the optimal choice for our case.

One of the most challenging tasks in SLAM is the egomotion, which contains the rotation and translation of the camera, estimation of the camera. Least Square Fitting is a non iterative registration algorithm firstly proposed by Arun et al. [START_REF] Somani Arun | Least-squares fitting of two 3-d point sets[END_REF]. The advantage of using Least Square Fitting algorithm is its nature of efficiency. However, since using Least Square Fitting algorithm consider all of the matching pairs, it is very sensitive to outliers. Another technique to estimate rigid motion between two point clouds is the minimum 3-Point RANSAC [START_REF] Jiang | Static-map and dynamic object reconstruction in outdoor scenes using 3d motion segmentation[END_REF] algorithm. This registration algorithm takes iteration to maximize number of inliers under the RANSAC framework [START_REF] Martin | Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography[END_REF]. The advantage of using 3-Point RANSAC is its nature to be robust to outliers.

Since the introduction of ICP by Chen and Medioni [START_REF] Chen | Object modelling by registration of multiple range images[END_REF] and Besl and McKay [START_REF] Paul | Method for registration of 3-d shapes[END_REF], the algorithm was further developed into various type of ICP by many researchers to make it more robust. Rusinkiewicz and Levoy [START_REF] Rusinkiewicz | Efficient variants of the icp algorithm[END_REF] nicely classified these variance of ICP into six main stages.

III. MATHEMATICAL MODEL

A. 3-Point RANSAC Point Cloud Registration

3-Point RANSAC algorithm can robustly and precisely estimate the rigid motion parameter R and t between two 3D point clouds in a linear manner. Let X = [x, y, z] T and Y = [x , y , z ] T be two corresponding 3D points under rigid transformation, we have:

X = RY + t, (1)
where R is the 3 × 3 rotation matrix and t is the 3 × 1 translation matrix. Let g be the Gibbs representation of rotation matrix R, we have G = [g] × is a 3 × 3 skew-symmetric matrix where g = e tan θ 2 with e = [e x , e y , e z ] T is the Euler rotation axis and rotation angle θ. Applying the Cayley Trasformation, R can be represented as:

R = (I 3 + G) -1 (I 3 -G), (2)
where I 3 is a 3 × 3 identity matrix. Replacing Eq. 1 using Eq. 2 and multiplying (I 3 + G) on both sides, we have:

X -Y = -G(X + Y) + t, (3)
where t = (I 3 + G)t = [ tx , ty , tz ] T . Eq. 3 is a linear system such that the parameters can be estimated using a Linear Least Square approximation.

To robustly find the good estimation, a Random Sample Consensus (RANSAC) [START_REF] Martin | Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography[END_REF] algorithm is adopted to maximize the number of inliers. Let N be the number of matching pairs that have the RMSE less than threshold τ , so called inliers, we aim to maximize the number of inliers by randomly select 3 sample matching pairs as candidates to estimate the transformation matrix:

N = argmax T Ξ n i=1 X i -TY i 2 < τ , (4)
where T = R t 0 1 is the desired transformation matrix in homogeneous coordinate. n is the total number of matching pairs. Symbol • 2 denotes the l 2 -norm. Notation Ξ n counts the number of pairs satisfying the criteria that the pointpair distance is less than threshold τ .

B. ICP Point Cloud Registration

Different from 3-Point RANSAC algorithm, the ICP algorithm establishes the matching pairs by considering the nearest point from the new model to the reference. Let X = {X 1 , ..., X j } be the reference point cloud, and Y = {Y 1 , ..., Y k } be the new model, a point to point ICP matching pairs can be established by:

X j ↔ Y k := min{ X j -TY k 2 }, (5)
with j ∈ [1, ..., m], k ∈ [1, ..., K],
where symbol ↔ denotes the point to point matching property, and symbol := can be interpreted as defined by. m and K are the maximum index of the reference and new model point cloud, respectively. Eq. 5 builds the matching pairs according to the closest point distance from reference point cloud to the new model. Since the ICP has complete point-to-point matching between two point clouds, its performance is highly depending on the level of noise or number of inliers.

E icp = argmin T 1 m m j=1 X j -TY j 2 , ( 6
)
where E icp is the normalize root mean square error distance (RMSE) required to be minimized. T is the refined transformation matrix after ICP.

IV. METHODOLOGY A. Point Correspondences

To search for feature correspondences in 3D space, ideally, 3D feature descriptors [START_REF] Luís A Alexandre | 3d descriptors for object and category recognition: a comparative evaluation[END_REF] should be used. However, in practice, matching the 3D points using 3D descriptors is difficult to be maintained due to lack of robust 3D feature descriptors [START_REF] Jiang | Static-map and dynamic object reconstruction in outdoor scenes using 3d motion segmentation[END_REF], [START_REF] Danda | 2d-3d camera fusion for visual odometry in outdoor environments[END_REF]. In contrast, 2D image feature descriptors is much more robust and reliable in the case of 2D-2D matching. Therefore, to overcome this problem, we take the advantage of the 2D-3D camera setup to infer the 3D-3D feature correspondences based on the detected 2D-2D correspondences.

In details, the point correspondences are not established directly between 3D feature points, but we infer them indirectly by building the 2D feature matching pairs from their corresponding image frames. Once corresponding points between two images are obtained, we associate the 2D features with their corresponding 3D points by the 2D-3D correspondences from the calibrated 2D-3D setup. Figure 3 illustrates how the 3D-3D correspondences are established through the inferences of 2D-2D correspondences.

B. Feature Matching

As discussed about, the SIFT or other feature detectors can detect only sparse feature points, and those feature points can be easily lost tracked during a long sequence tracking. Moreover, since the moving objects in the scene are relatively small, there are very few features can be detected from the moving objects in most cases. On the contrary, Optical Flow feature matching can be used instead, which produces very dense feature matching between two consecutive frames.

Remind that image sequence (video) is involved, OF feature is the optimal choice for our system. The main advantage of optical flow feature matching is that every pixel can be consider as a feature point, such that we have very dense features. Accordingly, features from moving objects are dense enough to be tracked and registered.

C. Feature Tracking

In practice, the point density of the 2D image is much denser than the 3D point cloud. A problem raises that not all the image points has 3D correspondences. Fig. 4 explicitly shows the problem that sparse 3D points corresponds to dense 2D image points. To detect the key frames automatically, a MSE threshold is defined such that new keyframes are chosen whenever the MSE of a local sequence is greater than the threshold.

Valid Features Projected 3D Points

2) Initialization and ICP Refinement: Sparsity and poor geometric structure (see Fig. 1) make ICP is not effective to register point clouds in our case, it will easily trapped in local minima. Even though 3-Point RANSAC can register the point clouds nicely, registration through transformation propagation leads to an error accumulation problem. It is already reduced by using pipeline in Fig. 5. However, to achieve a high quality registration, a refinement step is necessary. Therefore, combining 3-Point RANSAC as initialization and ICP for refinement can perform a very effective registration (see Fig. 8). Fig. 6: Key-Frame based Registration Pipeline 3) Recognizing Part of Vehicle: ICP refinement doesn't guarantee correct registration in all case. It fails in registering two point clouds with small partial overlapping (see Fig. 2). In this case, the ICP tries to minimize the registration error by rotating and translating the model, which produces incorrect registration.

To avoid that situation, the registration should be performed only for set of points that share the same side of the object (Figure 7a). Normal vector is able to give that information. Thus, by using normal vector, we are able to obtain those set of point.

To differentiate the faces of rigid object, a plane fitting algorithm performed. A plane is defined as z = ax + by + c. By using 3 non-colinear points, we can define a plane. To robustnessly fit the plane, the RANSAC framework is applied. Figure 7a and 7b shows the fitted back-side planes and leftside planes, respectively. Getting the fitted planes, only the points close to the fitted planes will be considered to perform registration.

V. RESULT AND EVALUATIONS

To evaluate the performances of the proposed algorithm, experiments are conducted on realistic outdoor environment by using KITTI dataset [START_REF] Geiger | Are we ready for autonomous driving? the kitti vision benchmark suite[END_REF]. The descriptions of the dataset are summarized in Table I. The software is developed using a computer with Intel Quad Core i7-2640M, 2.80GHz, 7.8GB Memory. To evaluate the performances with and without registration refinement, Mean Square Error (MSE) is chosen to quantify the performance. 3D point cloud from first frame is considered as a global reference model. And the registration error is defined as MSE of corresponding 3D points from registered model to the global reference model.

Figure 8 shows that both ICP refinement with and without interpolation (red and blue plot) have registration failure in around frame 30-35. These frames is when small partial overlapping occurred (see Fig. 9). Also, a denser point cloud interpolation is applied to increase the stability of the ICP registration.Satisfactory result appear for refinement with the proposed algorithm with plane fitting (green plot). After dividing different parts of vehicle by considering the normal consensus of point cloud subset, registration performed accurately, even for small partially overlapping cases.

Figure 10 shows qualitative results of overall registration from three different viewports. 3D point clouds from multiple frames are registered very satisfactorily. Comparing to single scan point cloud by 3D laser scanner (see Fig. 1 and Fig. 2) the the registered 3D point cloud of vehicle have much denser points and richer geometry structure details.

VI. CONCLUSION AND FUTURE WORK

In extreme cases, Iterative Closest Point (ICP) or its variance is not able to perform registration correctly. It easily trapped in the local minima. Highly sparse 3D point cloud, small partial overlap, existence of noise and inliers, create a new challenge to perform a registration, especially a registration dedicated to long sequence. Feature point reduction and error accumulation is the main problem in long sequence 3D point cloud registration. To overcome all of the problem defined, we took advantages from a calibrated setup of 2D-3D camera and proposed a algorithm for long sequence registration. The algorithm conducted in realistic outdoor environment by using KITTI dataset [START_REF] Geiger | Are we ready for autonomous driving? the kitti vision benchmark suite[END_REF]. From experiments, the proposed algorithm produced very satisfying results both quantitatively and qualitatively. However, the algorithm is still not sufficiently generic based on the assumption of plane fitting is adequate. To robustify the registration, robust estimation techniques, such as M-Estimator, can be applied to the algorithm. Furthermore, 2D images can also be employed as additional information to build up the Absolute Pose Estimation problem to estimate camera poses to register the point cloud.

Fig. 1 :

 1 Fig. 1: Illustration of Poor 3D Geometric Structure Problem: left image, the red and blue point clouds are acquired from two different scans. The right image shows the failure registration result by using ICP algorithm due to the poor geometric structure of the point clouds.

Fig. 2 :

 2 Fig. 2: Illustration of Small Partially Overlapping Problem: left image, the red and blue point clouds are acquired from two different scans. The right image shows the failed registration result using ICP algorithm.

Fig. 3 :

 3 Fig. 3: Inferring 3D-3D correspondences from 2D-2D feature correspondences.

  Features and 3D Projections Association in the NextFrame.

Fig. 4 :Fig. 5 :

 45 Fig. 4: Closest Feature Association Problem. In this figure, the white dots in both sub-figures represents the projections from 3D points. The green dots in subfigure (a) are the valid features after forward-backward OF

  (a) Similar side.

  on Plane (b) Back side plane.

  on Plane (c) Left side plane.

Fig. 7 :

 7 Fig. 7: Different Part the Rigid Vehicle.

FramesFig. 8 :

 8 Fig. 8: Comparison of Different Registration Algorithms.

  (a) Image from Frame 30. (b) 3D Model.

Fig. 9 :

 9 Fig. 9: Poor 3D Geometric Structure Frame.

  (a) Back-left View. (b) Back-right View. (c) Top View.

Fig. 10 :

 10 Fig. 10: Registration Results using Our Algorithm.

TABLE I :

 I Dataset Profile

	Name of Dataset	Van Dataset
	Number of Frame	44
	Number of Side	3
	Transition Frame	Left-Back Right-Back	1-30 31-44
	Number of Keyframe Produced	11


Ahmad Zawawi Jamaluddin
email: ahmad.jamaluddin@u-bourgogne.fr

Osama Mazhar

Olivier Morel

Ralph Seulin

David Fofi

Design and calibration of an omni-RGB+D camera

Keywords: fisheye, dioptric, stereo vision, unified model

INTRODUCTION

The research of bio-inspiration is increasingly important along with the creation of autonomous robotics to replace the manned system [START_REF] Marefat | A hemispherical omni-directional bio inspired optical sensor Electrical Engineering (ICEE)[END_REF]. The fabrication of the wide field of view (FoV) camera, work together with computer vision technique and integrated with the robotic system makes this research more effective. Nowadays, the camera is not only used for visualisation, but also maneuvering and analyzing the current state of the robot [START_REF] Li | Full-view spherical image camera[END_REF]. This paper will present the design and calibration of a new sensor composed of two fisheyes and a stereo-rig allowing both a 360 • field-of-view and a directional depth-map of the scene in one-shot.

HYBRID CAMERA SYSTEM

Omni-directional vision system

The omnidirectional vision can be achieved by using a single camera with special lens (dioptric) or a single camera with reflected mirror (catadioptric) also by combining several identical cameras (polydioptric) [START_REF] Nayar | Catadioptric omnidirectional camera[END_REF]. Fig. 1 shows the types of omnidirectional vision camera. Dioptric or fisheye camera has a better performance because it does not produce a blind area compare to catadioptric camera. The advantages of using two cameras are easy for synchronisation, low bandwidth, low power consumption and low cost.

Stereo Vision -ZED Camera

The 3D information from a digital image can be estimated by using a stereo vision system. This depth information are estimated based on the two-view geometry technique to approximate the 3D structure of the scene. ZED is a stereo vision camera developed by Stereo Labs [START_REF] Zed | [END_REF]. It has already been pre-calibrated by the factory.

PROPOSED SYSTEM AND RESULTS

Unified Camera Model

The intrinsic parameters of two fisheyes are calculated before projecting the image on the unit sphere. Mei's projection model [START_REF] Mei | Single view point omnidirectional camera calibration from planar grids[END_REF] has been used as a reference and its provides procedures to put the image on the unified model. The projection model is enhanced by incoporating it with the tangential and radial distortions. Let us consider a 3D point χ = (X, Y, Z) T on the camera frame, and project it to the unit sphere with C m as a center of sphere.

χs = χ χ = (Xs, Ys, Zs) T . (1)
Then the point χ s mapped to the new reference frame with the new center C p (χs)F m → (χs)F p = (Xs, Ys, Zs) T .

(2)

Next, the point projected onto the normalised plane

mu = ( Xs Zs + ξ , Ys Zs + ξ , 1) T = Xs. (3)
The models of distortion (tangential and radial) are added in the projection model. xc = x1 + k1r 2 + k2r 4 + k5r 6 + 2k3xy + k4(r 2 + 2x 2 ), (4) yc = y1 + k1r 2 + k2r 4 + k5r 6 + 2k4xy + k3(r 2 + 2y 2 ), [START_REF] Mei | Single view point omnidirectional camera calibration from planar grids[END_REF] where:

r = x 2 + y 2 , ( 6
)
and the sum of distortion is

m d = mu + D(mu, V ), ( 7
)
where V is the coefficient of distortion.

V = (k1, k2, k3, k4, k5, k6),

and finally, the point m d is projected to the image plane using K, which is a generalized camera projection matrix.

The value f and η should be also generalized to the whole system (camera+lens).

p = Km d =   f1η f1ηα u0 0 f2η v0 0 0 1   m d , (9)
where the [f 1 , f 2 ] T is the focal length, (u 0 , v 0 ) is the principal point and α is the skew factor. Finally, by using the projection model, the point on the normalized camera plane can be lifted to the unit sphere by the following equation:

-1 (mu) =      ξ+ √ 1+(1-ξ 2 )(x 2 +y 2 x 2 +y 2 +1 x ξ+ √ 1+(1-ξ 2 )(x 2 +y 2 x 2 +y 2 +1 y ξ+ √ 1+(1-ξ 2 )(x 2 +y 2 x 2 +y 2 +1 -ξ      . ( 10

Omni-Vision Camera Rig Setup

Our proposed system is a combination of two fisheyes and stereo vision camera. The two fisheye cameras with FoV of 185 • are fixed back to back. Thus, they cover the whole scene at once. A depth camera, namely ZED Camera is placed at 90 • to the right fisheye camera, allowing the system to observe the frontal view at higher resolution, with the depth information for each pixel of the image. Fig. 3 shows the proposed omni- vision camera rig. The calibration of each camera for intrinsic parameters is done using Omnidirectional Calibration Toolbox as described in [START_REF] Mei | Single view point omnidirectional camera calibration from planar grids[END_REF]. The extrinsic parameters (rotation and translation) of right fisheye camera with respect to the left fisheye camera is computed using a straightforward approach to match overlapping image features. The rig is placed in a sufficiently open space to acquire images simultaneously from the two fisheye cameras. The features are matched manually with raw images with lens distortion and the selected points are projected onto three dimensional space using the camera intrinsic parameters obtained from the calibration results. Rigid Transformation is estimated from a set of three dimensional points from each camera, thus rotation and translation are obtained as follows: The results show that the orientation of the right fisheye camera rotate about y-axis of 180 • . Assuming that both cameras have the same optical center, the translation obtained is neglected and both images from left and right fisheye cameras are projected onto a unit sphere to obtain an Omnivision of the scene as shown in fig. 4. The extrinsic parameters of the depth camera are obtained using a method described in [START_REF] Lbraly13 | Calibration of non-overlapping cameras-application to vision-based robotics[END_REF]. The omni-vision system is considered as a single view Fig. 4. Fisheye camera images projected and aligned onto the unit sphere vision system which provides a large scale SLAM. It provides the dense 3D reconstruction and texture mapping. Since the system has a prey and predator-like vision features, the system also provides sufficient information for target detection and tracking.

T =    -1.

CONCLUSION

We have successfully combined the two hemispheric images to the unit sphere by matching points along the overlapping area. This research is in progress and we expect its great potential in a practical applications.

Fig. 1 .

 1 Fig. 1. The omnidirectional camera: types and characterization (From Youcef Mezouar: Omnidirectional vision: Application to visual servoing).

) ξ 1 Fig. 2 .

 12 Fig. 2. Mei projection model from 3D world to 2D image. This figure is the courtesy of Christopher Mei [5].

Fig. 3 .

 3 Fig. 3. Proposed omni-vision camera rig.

ACKNOWLEDGMENTS

This work was supported by the Regional Council of Burgundy (CRB) and FEDER (grant nb. 2015-9201AAO050S02760).