Katia Jaffrès-Runser

Qi Wang

Jean-Marie Gorce

Cristina Comaniciu

Fundamental limits of wireless ad hoc networks: upper MO bounds

Keywords: fundamental limits, capacity, delay, energy, wireless networks, upper bound, multiobjective bound limites fondamentales, capacité, délai, énergie, réseau sans-fil, borne supérieure, borne multicritères

This report addresses the problem of deriving fundamental trade-off bounds for wireless ad hoc networks when multiple performance criteria are of interest. It proposes a MultiObjective (MO) performance evaluation framework composed of a broadcast and interference-limited network model, a steady state performance metric derivation inspired by a discrete Markov chain formalism and formulates the associated MO optimization problem. Pareto optimal performance bounds between end-to-end delay and energy for a capacity-achieving network are given for the 1-relay and 2-relay networks and assessed through simulations.

Introduction

Wireless ad hoc or sensor networks are many times operating in difficult environments and require several performance criteria to be satisfied, related to timely, reliable, and secure information transfer. Routing and resource allocation protocols are key solutions for ensuring proper information transfer across the network but their design has to cope with imperfect wireless transmissions. Finding the fundamental performance limits of wireless ad hoc networks is a very active research area (cf. [2] and references herein) that has triggered a comprehensive research effort in the last decade. Initial research has focused on defining the capacity region of the network [3,[START_REF] Toumpis | Capacity regions for wireless ad hoc networks[END_REF] and finding capacity maximizing resource allocation and routes [5].

Depending on the application, additional metrics such as delay, outage or energy consumption have become relevant. They are inter-dependent and can not be optimized separately. In this perspective, new multiobjective frameworks (network models and optimization tools) are needed to derive a multiobjective (MO) performance bound defined as the fundamental trade-off between the objectives of interest. These MO bounds could provide insight to improve network design and performance, as well as an upper bound against which to compare the performance of existing protocols [2].

The derivation of bounds relies on the definition of cross-layer network models capturing wireless variations on which a joint optimization of routing and resource allocation is performed [START_REF] Toumpis | Capacity regions for wireless ad hoc networks[END_REF]5]. As proposed by Goldsmith et al. in [2], a promising way towards characterizing fundamental MO bounds in wireless ad hoc networks is to exploit "the broadcast features of wireless transmissions through generalized network coding, including cooperation and relaying". In this paper we propose a framework composed of a cross-layer network model and a steady-state performance evaluation model capturing capacity, delay and energy metrics. We formulate the associated MO optimization problem whose resolution provides both the MO bound and the MO optimal network configurations. Unique to this work is that our network model and MO problem has been designed to incorporate broadcast and interference-limited channels and thus, is capable of deriving bounds for a layerless communication paradigm [2] that integrates generalized network coding, cooperation and relaying.

The paper is organized as follows. The descriptions of our cross-layer network model and the relative MO optimization problem are given in Section 2 and 3, respectively. Section 4 presents our steady-state performance evaluation metrics. First results are presented in Section 5 and Section 6 concludes the paper.

Network model

We assume a synchronized network where transmissions are time-multiplexed. A frame of |T | time slots is repeated indefinitely. One or more packets (or symbols) can be transmitted in a time slot. In the rest of this paper, our examples assume that one packet is being sent in one time slot. A time epoch s is defined as the time needed to transmit one frame of |T | time slots. For any time slot u ∈ T , there is an interference-limited channel between any two nodes i and j of the network. This channel is modeled by the probability of a symbol or a packet to be correctly transmitted between i and j in time slot u. This probability is referred to as the channel probability and denoted p u ij in the following. It models interference as an additive noise and is computed considering the distribution of the bit error rates (BER) or the packet error rates (PER) as shown hereafter. In this graph, an edge (i, j, u) represents the channel between nodes i and j in time slot u. Each edge is assigned a weight of p u ij . For each node i ∈ V, -→ N i and ← -N i are the set of edges leaving from and going into i, respectively. Each channel is assumed to be in a half-duplex mode, i.e. a node cannot transmit and receive a packet at the same time. An emission is defined as the couple (i, u) ∈ V × T and represents the fact that node i is emitting in a time slot u.

X i p 1 ij p u ij τ u i τ 1 i τ |T | i Node j Node i p |T | ij X j
A set of sources O and destinations D is defined. A source S o ∈ O, o ∈ {1, .., |O|} can emit a flow of data to a single destination D d ∈ D, d ∈ {1, .., |D|} or to several ones. We make the assumption that source and destination nodes do not relay the information. Sources only emit packets and destinations only receive packets. Multi-hop transmissions are allowed and we model the other nodes as relay nodes R = V -O -D. We have N = |R| the number of relays in the network.

The network is synchronized. Depending on their time slot assignments, sources and relays emit their symbols (or packets) at the beginning of their assigned time slots. Nodes emitting in the same time slot may interfere, reducing reception probability of their emitted data. Nodes that are not emitting in a time slot can receive symbols (or packets) in this time slot. We assume that relays have an incoming buffer and |T | outgoing buffers. All buffers are able to store the amount of symbols (or packets) transmitted in one time slot duration.

In our examples, they can store one packet. A relay receiving a packet has to decide if it will discard it or in which slot of the next frame it will send it. We consider as well in our model that a relay can not differentiate packets: identical packets are indiscernible.

Emission rate

A node j emits a flow of symbols (or packets) in time slot v. Symbols are the realization of a random variable -→ X v j chosen from an alphabet X (X = {0, 1} for instance). In the graph, it means that node j emits the same flow on its Inria -→ N v j outgoing edges using the time slot v. Thus, our model incorporates the broadcast property of wireless communications.

A flow of symbols coming into node j from node i on time slot u is modeled as a random variable ← -Y u ij . The symbols are transmitted by node i on the edge (i, j) ∈ E in time slot u and consequently, experience the probability p u ij of being received successfully in j. Depending on the level of interference, symbols are received successfully or not at j.

Having this, let ← -

Y j = { ← - Y ij , (i, j) ∈ ← - N j }
be the random variable giving all the symbols that can be received on all incoming edges of node j where ← -

Y ij = { ← - Y u ij , (i, j) ∈ ←- N u j }. Let -→ X j = { -→ X v j , j ∈ -→ N v j }
be the random variable giving all the symbols that can be emitted on the outgoing channels of node j for all the time slots. The relation between the outgoing random variables -→ X j and the incoming random variables ← -Y j defines a coding scheme for the network. In general, a relay has one main decision to take upon receiving a symbol or a packet: whether it should forward it or not. If it decides to forward it, the next step is to decide on which time slot to transmit it. If network coding is considered, the choice of the code adds to the set of decisions.

These decisions strongly influence the quality of the transmission and aim at mitigating the transmission errors due to fading and interference. Depending on these decisions, the emission rate of a node i on a time slot u varies. For instance, if a node decides to drop one packet out of two received on a same time slot u, the emission rate on channel u would become half the rate at which it received packets on the same time slot. Similarly, it is possible that a node transmits one packet every two received packets on time slot u because it is applying a coding scheme for which two packets are combined into a single transmitted packet through network coding. In both cases, the node is adapting the emission rate on each channel to combat errors.

Based on this statement, we characterize the behavior of a node i by the rate at which it is emitting its flow in each time slot u. The emission rate of node i in time slot u is denoted by τ u i . In our model, it is normalized and thus its value belongs to the interval [0, 1]. The emission rate of destinations is equal to zero. Having this, a vector of emission rates for each time slot can be defined

τ i = τ 1 i . . . τ |T | i
. Emission rate τ u i can be interpreted as the probability node i is emitting packets in time slot u. Assume a time slot permits the transmission of one packet, a node i transmitting at a rate of τ u i = 0.5 in time slot u will decide with probability 0.5 to transmit its packet in time slot u in the upcoming frame. If u∈T τ i = 0, node i isn't active.

The average rate at which all the symbols are coming into node j is given by ← -

r j = u∈T (i,j)∈ --→ N u j τ u i p u ij for p u ij
the probability for the symbol to be correctly received. The rate at which symbols are being transmitted by node j is -

→ r j = v∈T τ v j . Let τ = τ † 1 • • • τ †
N † be the emission rate matrix. A particular instance of τ belongs to the set of possible emission rate matrices denoted as Γ. An instance τ ∈ Γ is feasible if Properties 1 and 2 hold for each node:

Property 1: Flow conservation. The sum rate of all outgoing flows is lower or equal to the rate of all incoming flows, i.e.

- → r j ≤ ← - r j , ∀j ∈ V (1)
Property 2: Half duplex. A node j is able to receive a message on a time slot u the proportion of time it is not transmitting on that same time slot. As a consequence, a τ is feasible if there is enough time to transmit and receive in any time slot for each node:

← - r u j + τ u j ≤ 1, ∀(j, u) ∈ V × T (2)
where ← -

r u j = (i,j)∈ --→ N u j τ u i p u ij stands
for the incoming cumulative rate in time slot u.

We define a set of active emissions A corresponding to the values of τ . An emission (i, u) ∈ V × T is said to be active if τ u i > 0. As a consequence, the set A is given by

A = {(i, u) ∈ V × T | τ u i > 0}.
Similarly, the set A u refers to the set of active emissions restricted to time slot u. We have

A u = {(i, v) ∈ V × T |τ u i > 0, v = u}

Channel probability

Knowing τ , it is possible to derive the interference level for each edge of the network. We want to stress here that the channel probability derivation is computed for a steady state where all nodes with τ u i = 0 are transmitting concurrently in time slot u. We model interference as an additive noise computed at the end of each edge (i, j, u) ∈ E of the network. We recall that p u ij is the probability for a symbol emitted by node i to arrive successfully at node j in time slot u. It is a function of the statistical distribution of the Signal to Noise and Interference Ratio (SINR) at the location of the destination node j. This quantity can either be defined using a Packet Error Rate (PER) or a Bit Error Rate (BER) depending if transmissions on the network are packetized or not. In the following we consider that transmissions are packetized. Before giving the exact expression of the channel probability, a few preliminary definitions and notations are given hereafter:

Pathloss attenuation factor and transmission power a ij reflects the attenuation due to propagation effects between nodes i and j for a simple isotropic model. All nodes use the same power P T .

Interference Since transmissions are time-multiplexed, interference only occurs between transmissions using the same time slot. Let I u ij be the power of interference on the edge (i, j, u) ∈ E computed at node j. If I u ij is a set of nodes interfering the emission of i at node j in time slot u, interference power is calculated as

I u ij = k∈I u ij P T • a kj .
SINR SINR between nodes i and j in time slot u is defined by

γ u ij = P T •aij N0+I u ij
where I u ij is the interference power on the link and N 0 the noise power density.

Packet error rate (PER) For a specific value of SINR γ, the packet error rate P ER is computed with

P ER(γ) = 1 -[1 -BER(γ)] N b
where N b is the number of bits of a data packet and BER(γ) is the bit error rate for the specified SINR per bit γ which depends on the physical layer technology and the statistics Inria of the channel. Results are given for an AWGN channel and a BPSK modulation without coding where BER(γ) = 0.5 * erfc( √ γ).

We assume that there is no medium access mechanism and that emissions are totally independent. Within a time slot u, channel access is purely randomized and depends on the emission probabilities of the active nodes. Thus, on edge (i, j, u), the set of interferers present is random. In order to derive an average channel probability over time, we compute the average power of interference on the edge over time. Therefore we have to enumerate all possible configurations where active nodes interfere with the emitter i. These configurations are designated as interfering sets.

Interfering sets We recall that A u is the set of active emissions in time slot u. An interfering set I u ij for edge (i, j, u) belongs to the set of all possible interfering sets L u ij on this edge. L u ij is the power set P(A u -{i}), i.e. the set of all subsets of A u -{i}. Because of the half-duplex constraint, the receiver j is kept in the set A u -{i}. Thus, the channel probability computation accounts for the interfering sets where j is active. If j is active, the SINR is very low and transmission on the edge (i, j, u) is impossible. Equation ( 3) details the derivation of the channel probability p u ij as the average of the PER experienced for all possible interfering sets l ∈ L u ij on edge (i, j, u) referred to as P ER l :

p u ij = l∈L u ij [1 -P ER l ] .P l (3)
where P ER l = P ER(γ l ), l ∈ L u ij and γ l is the SINR experienced on the edge (i, j, u) when the nodes of the interfering set l are active.

P l is the probability for nodes of the interfering set l to be active and create interference on (i, j, u). More specifically, it is the probability that the nodes of the interfering set l are transmitting concurrently and the others are not:

P l = k∈l τ u k • m∈{A u \l} (1 -τ u m ) (4)
In ( 4), k∈l τ u k gives the probability that the active nodes of the interfering set l are transmitting and the other product the probability that the {|A u \ l}| other active nodes are not.

Similarly to the vector of emission rates, we define the vector of channel probabilities for a link (i, j) ∈ E:

P ij = p 1 ij • • • p |T | ij . Let P j = P † 1j • • • P † N j †
be the matrix giving all incoming channel probabilities at node j and P = P 1 • • • P N the matrix of all channel probabilities.

Forwarding and scheduling decisions

This section introduces forwarding and scheduling decisions of the nodes. These decisions are represented by the probability x uv ij of a node j to transmit on time slot v a packet coming from node i on time slot u. We will refer in the rest of the text to the forwarding probability x uv ij . For each node of the network, we can define a N |T |-by-|T | matrix giving all the forwarding probabilities relative to any node j of the network as follows. This forwarding matrix is given by

X j = X 1j • • • X N j †
where each matrix X ij provides the scheduling probabilities of a flow of packets coming from node i on its output times slots, depending on the time slot the packets are received on. We have

X ij =     x 11 ij • • • x 1|T | ij . . . . . . x |T |1 ij • • • x |T ||T | ij     .
The matrix of forwarding probabilities is related to the matrix of emission rates τ and the matrix of channel probabilities P with the following set of |A| equations

(i,j)∈ --→ N u j u∈T τ u i p u ij x uv ij = τ v j , ∀(j, v) ∈ A (5)
where τ u i p u ij is the probability that a packet sent by i on time slot u arrives in j. These equations ensure a flow conservation. They strictly constrain the choices of forwarding probabilities.

The forwarding probabilities represent the decisions of the nodes to either (i) retransmit all the packets or symbols received or (ii) reduce the output rate by dropping or re-encoding them together. From now on, we will refer to the set of all forwarding probabilities of the complete network using a matrix

X = [X 1 . . . X N ] , X ∈ X of size N.|T |-by-N.
|T | where X is the set of all possible matrix instances.

Illustrative examples

1-relay network In this example, we consider a basic 1-relay network as the one depicted in Fig. 2 (N = 1). Location of all nodes are known. For clarity S 0 and τ R = 0 τ 2 R . The forwarding probabilities matrix is:

D S p SD p RD τ S X R τ R p SR
X R = x 11 SR x 12 SR x 21 SR x 22 SR = 0 τ 2 R /(τ 1 S .p 1 SR ) 0 0
where the values of the x uv ij are derived according to the system of ( 5). From our model, we can for instance easily derive the average number of packets arriving at D which is the sum of the probabilities to receive a packet from all possible paths: AB , x 32 BA , x 12 SA et x 13 SB . The system of ( 5) provides 2 constraints on the forwarding probabilities:

f = τ 1 S .p 1 SD + τ 1 S .p 1 SR .x 12 SR .p 2
τ 1 S .p 1 SA .x 12 SA + τ 3 B .p 3 BA .x 32 BA = τ 2 A τ 1 S .p 1 SB .x 13 SB + τ 2 A .p 2 AB .x 23 AB = τ 3 B (6)
For a set of forwarding probabilities that respect the constraints in [START_REF] Toumpis | Capacity regions for wireless ad hoc networks[END_REF], the average number of copies f of a packet emitted by S arriving at D is computed by the infinite sum of the probabilities of receiving a packet through any possible path. In the two relay network, there is a loop between A and B which may create an infinite number of paths of increasing length. If p 2 AB x 23 AB p 3 BA x 32 BA < 1, f can be computed by:

f = τ 1 S p 1 SD + 1 1-p 2 AB x 23 AB p 3 BA x 32 BA p 1 SA x 12 SA p 2 AD + p 1 SB x 13 SB p 3 BD + p 1 SA x 12 SA p 2 AB x 23 AB p 3 BD + p 1 SB x 13 SB p 3 BA x 32 BA p 2 AD

MO optimization problem

The scope of this section is now to take advantage of the previously described network model to derive a framework capable of extracting the set of Paretooptimal networking solutions with respect to given performance criteria. A Pareto optimal set is composed of all the non-dominated solutions of the MO problem with respect to the performance metrics considered. The definition of dominance is: Definition 1: A solution x dominates a solution y for a n-objective MO problem if X is at least as good as y for all the objectives and x is strictly better than y for at least one objective.

Solution and search space

Based on our network model, several parameters can be treated as optimization variables: the location of the N relays, the number of relays N , the number of time slots |T | and the forwarding probabilities represented by matrix X ∈ X . In this paper, considered variables are the location of the N relays and their respective forwarding probabilities. Location of relays can be chosen in a convex set C. For one network realization, the forwarding probabilities are chosen in the set X . Thus the complete search space is S = C N × X .

Let f C , f D and f E be the performance criteria relative to capacity, delay and energy, respectively. Capacity is maximized, while energy and delay are minimized. The derivation of these metrics is provided in the next section. Our goal is to solve the following multiobjective optimization problem by finding the set of Pareto-optimal solutions S opt :

S opt = {x ∈ S | ∀x ∈ S opt , ∀y ∈ S c opt , x y} (7)
where S opt ∪ S c opt = S and notation x y means that x strictly dominates y. The cardinality of the search space S as defined grows exponentially with N ,|T | and |C|. There is a whole set of constraints for a solution x ∈ S to be valid already presented in (5). As a consequence, the solutions of S that do not follow these constraints are dropped before being evaluated by the search algorithm.

The question of choosing the optimal value N is out of the scope of this paper but will constitute an important extension of this work. This value is a trade-off between the search complexity and the quality of the Pareto-set obtained. Here, we consider small values for N . Similarly, results are presented for small values of |T | since we can achieve interference free transmission for |T | = N + 1 time slot allocations.

Change of variable

The derivation of the transition rate matrix τ ∈ Γ knowing a solution X ∈ X is intractable. In a nutshell, to compute the emission rates of a node j, X and the incoming transmission probabilities P j have to be known yet the elements of P j are a function of the emission rates matrix, creating a circular dependency between the variables.

Nevertheless, being able to model exactly the interference level of the network is appealing and to keep the same formalism, we propose a reverse approach where a solution of the optimization problem is defined by the set τ of the emission rates for all the nodes and time slots. From τ , it is possible to derive the channel probabilities matrix P according to (4) since the activity of all nodes on each time slot is known.

Only instances of τ that meet the constraints relative to Property 1 and 2 are further considered as valid. Now that we have a valid τ , we can derive all the forwarding matrices X ∈ X that verify the constraints of equation (5). There are |A| constraints, each one constraining the choice of the x uv ij for all nodes and time slots of the network with respect to τ . Let X τ be the subset of X that verifies (5) with respect to the emission rate matrix τ . Each solution X ∈ X τ can be evaluated according to f C , f D and f E . The MO optimization problem of (7) stays unchanged, however, the way the search space is constructed for one network configuration has changed and permits to select in a first stage a subset of valid solutions X τ . This section provides a framework for the definition of performance criteria for a particular solution x ∈ S of our MO optimization problem. It is an important contribution because it permits a fast performance derivation for end-to-end criteria such as capacity, delay and energy. The proposed derivation relies on the definition of a transition matrix which is composed of the probabilities for a flow of packets to be re-transmitted by the nodes of the network. Such a formulation is inspired by the theory of Markov chains. However, we want to stress that we do not model the network using a Markov chain. We define instead a transition matrix that has the properties needed to be able to re-use some results from the theory of Markov chains. Table 1 gives notations specific to this section.

Inria

Flow vectors

At time epoch 0, a source S o ∈ O emits its flow of packets in a frame of |T | time slots following its emission rate vector. Here we make the assumption that the network is already fully loaded and that all relays and sources are active following the emission rate matrix τ . After t time epochs, the flow of S o has travelled in the network. All or part of it has reached different relays and because of losses on the channel, routing and scheduling decisions, its rate has most probably evolved. Let -→ F u i (s) represent the rate at which this flow is being emitted by node i in time slot u after s time epochs. Similarly, ← -F u i (s) represents the rate at which the same flow is being received at node i in time slot u from all its neighbors after s time epochs. For the rest of the paper, we will loosely designate -→ F u i (s) and ← -F u i (s) as outgoing and incoming flow rates respectively. The incoming flow rates are mostly of interest for destinations.

F (s) represents the vector of outgoing flow rates at the relays and the incoming flow rates of the destinations for each time slot after s time epochs. F (s) is referred to as the flow rate vector and is given by:

F (s) = -→ F R (s) ← - F D (s)
where -→ F R (s) stands for the outgoing flow rates emitted by the relays and ← -F D (s) for the incoming flow rates received at the destinations after s time epochs. More specifically, we have: The flow rate vector is of dimension (N +|D|)|T | if all possible transmissions are accounted for. However, if node i never transmits in a time slot u (i.e.

-→ F R (s) = -→ F 1 (s) • • • -→ F N (s) with -→ F i (s) = -→ F 1 i (s) • • • --→ F |T | i (s) and ← - F D (s) = ← -- F D1 (s) • • • ← --- F D |D| (s) with ← -- F D d (s) = ← -- F 1 D d (s) • • • ←-- F |T | D d (
τ u i = 0), -→ F u i (s) = 0.
As a consequence, the flow rate vector can be reduced to the size of the active emissions. Hence, the size of F (s) can be reduced to

(|A| + |D|)|T |.

Transition matrix

The propagation of the packets on the edges and the forwarding decisions of the nodes can be represented by a transition matrix M . When it is applied to the flow rate vector of time epoch s, the new flow rate vector of time epoch (s + 1) is obtained. Formally we have:

F (s + 1) = F (s) • M
The transition matrix M has the following structure: Transition matrix M has a similar canonical structure than a finite absorbing Markov chain [4]. Here, we have l transient states, i.e. relay nodes that are forwarding packets using active transmissions. The relaying matrix Q gives the probabilities for any emission (i, u) at time epoch s to be emitted as (j, v) at the following time epoch (s + 1) by the relays of the networks. We have m absorbing states given by the reception of packets at the destinations in D. Identity matrix represents the fact that packets received by destinations are never forwarded, but absorbed. Matrix D is composed of the probabilities to go from any transient state to any absorbing state, i.e. the probabilities for any emission (i, u) at time epoch s to arrive at a destination D d in time slot v at time epoch (s + 1) .

The relaying matrix Q is structured as follows:

Q =      0 Q 12 • • • Q 1N Q 21 0 • • • Q 2N . . . . . . Q N 1 • • • Q N -1N 0      0 is an |T |-
Q ij =     Q 11 ij • • • Q 1|T | ij . . . . . . Q |T |1 ij • • • Q |T ||T | ij     (8)
where Q uv ij is the probability for a node j to retransmit on channel v a packet that has been transmitted by node i on time slot u. From our network model, it equals:

Q uv ij = p u ij x uv ij
The arrival matrix D is given by: The channel probability computation of (3) accounts for the interference distribution created by all sources and relays in the network. Thus, once the channel probabilities are set, it is possible to derive the outgoing flows for the packets of each individual source using the same source and transition matrices. If follows that M S and Q can be applied to any type of source-destination communication (unicast, multicast, ...). This is a direct consequence of the fact that the inherent broadcast and interference-limited properties of the radio channel are accounted for in our network model. For instance, for a multicast transmission, the arrival matrix D is defined for |D| > 1 destinations. The matrix Q stays unchanged compared to a unicast transmission with |D| = 1. In a convergecast scenario where multiple sources are transmitting to a same destination (a sink for instance in a wireless sensor network), each source is propagated in the network using the same Q to a unique destination modeled by D. An important consequence for our optimization problem, is that its complexity is not a function of the number of flows transiting in the network.

D =    D 1D1 • • • D 1D |D| . . . . . . D N D1 • • • D N D |D|    D iD d is a |T |-
In the following, matrices are represented considering all N relays. However, when evaluating a specific solution τ , these matrices are derived only for the set A of active emissions and the size of the square matrix Q is reduced from l = N.|T | to l = |A|.

Fundamental matrix

We really want to stress that the proposed transition matrix does not follow the definition of a regular Markov chain. Indeed, the sum of the probabilities on a line of M is greater than 1. However, if Q is convergent, Q s tends to zero as s tends to infinity and the transmission reaches a steady state. In this case, the following theorem of absorbing Markov chains still holds:

Theorem 1: I -Q has an inverse, and

(I -Q) -1 = I + Q + Q 2 + • • • = ∞ s=0 Q s
Proof A complete proof is omitted for conciseness purposes but main intuition is given. For Theorem 1 to hold, Q must be convergent. In absorbing Markov chain theory, this is the case if from any non-absorbing state it is possible to reach an absorbing state. An element of Q gives the probability of a node to re-emit a packet transmitted by another node. Due to the broadcast property, a packet can be re-emitted by multiple nodes, introducing non-null elements in the relay flow vector -→ F R (s) as the number of time epochs (and hops) increases. Thus, any created packet (original or copy) has to reach a destination node to be absorbed. This is true for any packet traveling on a path without loop since it is being absorbed by a destination.

Packets traveling on paths with loops may infinitely circle in the loop and never reach a destination. We can prove that for the elementary loop present in the 2-relay network, packets are absorbed if the maximum eigenvalue of Q is strictly lower than one. This is the case if

Q 23 AB •Q 32 BA < 1, i.e. p 2 AB x 23 AB p 3 BA x 32
BA < 1 which is the probability for a packet entering the loop at A to be received at A again. For a larger loop, it can be shown by induction that if the probability of a packet entering the loop at a node i to be received at i (which we designate as the loop probability) is higher than or equal to one, Q is not convergent. Thus, Q is convergent if all paths with loops have a loop probability strictly lower than one.

Practically, solutions x ∈ S with non convergent Q can be disregarded from the MO optimization search since they lead to solutions creating an infinite number of received packets, which is detrimental to energy and delay minimization.

The inverse of I -Q is defined as the fundamental matrix and denoted M F = (I -Q) -1 . An element M F uv ij of the fundamental matrix represents the average number of re-emissions node j has done in v from packets sent by i in u until all packets are received at the destinations.

Knowing M F , it is possible to derive from F (1) the steady state flow vector

F (∞) = -→ F R (∞), ← - F D (∞) : -→ F R (∞) = -→ F R (1) • M F Inria ← - F D (∞) = -→ F R (∞) ← - F D (

Optimization criteria

The definition of our optimization criteria are directly derived from the fundamental matrix. Indeed, if I -Q is invertible, M F gives a measure of the performance of one solution x ∈ D when the number of time epochs tends to infinity. The very interesting feature is that it also accounts for all the possible cycles in the graph if Q is convergent and consequently models the broadcast property of the wireless channel. The criteria here after are defined for one source-destination flow. In the following, we consider source S o ∈ O transmits to destination D d ∈ D. Since the underlying network model holds for concurrent flows, it is possible to easily define criteria for multiple flows in the future.

Redundancy, capacity and reliability The number of packets absorbed at the destinations for one packet emitted by S o is given by:

f = ← - F D (∞) • I D (d).
I D (d) is a selection vector of dimension |D||T | where the |T | elements relative to destination D d are equal to 1 and the others are equal to 0. This vector accumulates the packets received in each time slot at destination D d . This criterion means that, with network solution x ∈ D, f packets are received at D d for one packet sent by S o . However, since relays are not able to discriminate packets, the set of f packets may be composed in the worst case of f copies of a packet sent by the source. What we can say here is that:

• if f ≥ 1,
there may be at most f different packets. Consequently, f provides an upper bound on both network capacity and packet redundancy. Let f C be the bound on the real capacity of the network:

f C = min(1, f )
This bound is exact if there exists a practical transmission scheme where all packets received at the destination are different.

• if f < 1, f can be interpreted as a reliability criterion defined as the probability to get one packet at D d .

Delay criterion Let P (H = h) be the probability for a transmission towards D d to be done in h hops. After s time epochs, a packet can travel up to h = s+1 hops. Thus

P (H = h) = ← - F D (1) • I D (d) h = 1 -→ F R (1) • Q h-2 • D • I D (d) h ≥ 2
We assume that a relay introduces a delay of 1 unit. Consequently, a h-hop transmission introduces a delay of h -1 units. The average end to end delay is computed by the infinite sum:

f D = ∞ h=1 (h -1) • P (H = h) = 0 • ← - F D (1)+ -→ F R (1) • I + 2Q + • • +(h + 1)Q h + •• • D • I D (d)
Let V be equal to

(I + 2Q + • • +(h + 1)Q h + ••). We can show that V = (M F ) 2 and consequently f D = -→ F R (1) • (M F ) 2 • D • I D (d) Proof Infinite sum ∞ i=0 i • Q i is equal to Q (1-Q) 2 . This sum is equal to Q • V . Thus, V = 1 (1-Q) 2 . By definition of M F = (1 -Q) -1 , V = (M F ) 2 .

Energy criterion

The energy is measured by the number of packets being emitted in the network by the source and the relays. By definition of the fundamental matrix, it is

f E = |T | t=1 τ o + -→ F R (1) • M F • 1 • • • 1 †
where 1 • • • 1 is a vector of N |T | ones that accumulates all the packets sent by the relays that have participated in the transmission. The first sum of the equation counts the number of packets sent by the source S o .

Illustrative examples

1-Relay network We still consider the example of Fig. 2. The initial flow vector F (1) is decomposed as

-→ F R (1) = 0 τ 1 S .Q 12 SR and ← - F D (1) = τ 1 S .D 1 SD 0 . Since N = 1, M F = I and -→ F R (∞) = -→ F R (1)
.

By definition, f = -→ F R (∞). 1 1 = τ 1 S p 1 SD + τ 1 S p 1 SR x 12 SR p 2 RD . Criterion f D = τ 1 S .p 1 SR .x 12 SR .p 2 RD since f D = -→ F R (∞) • 0 0 0 D 2 RD • 1 1
and

f E = τ 1 S + τ 1 S .p 1 SR .x 12 SR since f E = τ 1 S + -→ F R (∞) • 1
-→ F R (1) = τ 1 S Q 12 SA Q 13 SB and ← - F D (1) = τ S D 1 SD 0 0 . Q and M F are equal to: Q = 0 Q 23 AB Q 32 BA 0 ; M F = 1 |I -Q| 1 Q 23 AB Q 32 BA 1 with |I -Q| = 1 -Q 32 BA .Q 23 AB .
Steady state outgoing relay flows are:

-→ F R (∞) = τ 1 S |I -Q| Q 12 SA + Q 13 SB .Q 32 BA Q 13 SB + Q 12 SA .Q 23 AB and the arrival matrix is D = 0 D 2 AD 0 0 0 D 3 BD .
The three criteria are:

f = G + τ 1 S |I -Q| .(E + F ) with E = (Q 12 SA +Q 13 SB .Q 32 BA ).D 2 AD , G = τ S D 1 SD and F = (Q 13 SB +Q 12 SA .Q 23 AB ).D 3 BD . f D = τ 1 S |I -Q| 2 (A + B) with A = D 2 AD [Q 12 SA (1 + Q 23 AB .Q 32 BA ) + 2Q 13 SB .Q 32 BA ] and B = D 3 BD [Q 13 SB (1 + Q 23 AB .Q 32 BA ) + 2Q 12 SA .Q 23 AB ]
.

f E = τ 1 S + τ 1 S |I -Q| (Q 12 SA + Q 13 SB .Q 32 BA + Q 12 SA .Q 23 AB + Q 13 SB ).

Results

In this section, we look for the Pareto bounds that concurrently minimize capacity-achieving delay f * D and energy f * E criteria for the 1-relay and 2-relay problems. Here, capacity achieving criteria f * D and f * E are obtained by dividing the value of f D and f E by min(f C , 1) respectively. It permits to incorporate the effect of limited capacity on delay and energy in our analysis. For example, if f = 0.5, f * D = 2f D and f * E = 2f E which means that 2 times more packets have to be sent in average to reach perfect capacity which incurs double delay and energy. Theoretical Pareto-optimal solutions and bounds are obtained using non-dominated sorting genetic algorithm (NSGA-2) [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF]. To assess our network model, we simulate the Pareto-optimal solutions using the event-driven network simulator WSNet 1 .

Simulation settings

The distance between S and D is 620 meters. For P T = 0.151mW and a pathloss exponent of 3, direct transmission probability without interference is For NSGA-2, we use a population of size 500 and 1000 generations. The crossover probability is set to 0.9. Each Pareto optimal solution of the MO bound is simulated with WSNet. In our simulations, a perfect TDMA is implemented following 1-relay and 2relay problem specifications. S sends a packet every first time slot of every frame. Experiment is run for 10000 frames. For each Pareto optimal solution, the location of the relays and their forwarding probabilities are known. These forwarding probabilities are used in simulation as the decision to broadcast a packet upon its reception.

Values of f * D and f * E are calculated from simulated values of f , f D and f E derived as follows. f is measured by the total number of packets N rx received at D (including copies) divided by the number of packets transmitted by S. f D is the average transmission duration of the packets received at D. It is measured using the statistical distribution of the delays of the packets arrived over each possible distance measured in hops: P (h) = n(h)/N rx , where n(h) is the number of packets arrived in h hops at D. f D is then calculated with f D = ∞ h=1 (h -1).P (h). In this calculation, we consider as in our model that a transmission on a path of h hops introduces a delay of h -1 units. f E is the sum of the number of packets transmitted by the source and the relays.

Error between our model and simulations is measured for each criterion with a normalized root-mean-square error (RMSE) computed as

RM SE = 1 N opt N i=1 (f model (i) -f sim (i)) 2 f model (i) 2
where f model and f sim are the model and simulated values for criterion f and N opt is the total number of Pareto-optimal solutions in the MO bound.

Simulation results

Table 2 gives the values of RMSE for both 1-relay and 2-relay problems. Values are really small, showing a quasi-perfect match between the model and simulations. Fig. 4 plots the delay-energy Pareto bounds for both 1-relay and 2-relay problems. Delay results of 1-relay and 2-relay cases have been divided by |T | = 2 and |T | = 3 respectively to be comparable. A clear compromise is visible: decreased energy is obtained at the price of increased delay. Less reliable solutions spend less energy but need more transmissions to achieve perfect capacity.

For the 1-relay case, Pareto-optimal relays are located in a small circle of radius 10 -1 m located around the middle of [S, D] and have a forwarding probability x 12 SR = 1. Variability of optimal criteria is small ( 8e10 -2 ), meaning that the most important variable is here the location of R. The 2-relay bound dominates the 1-relay bound, thus transmission is more efficient using 2-relay paths than 1-relay paths. It is a consequence of the path diversity introduced by the 2 relays. Solutions with minimal delay have relays close to each other in the center of [S, D] and no loop exists (x 23 AB = 0 OR x 23 AB = 0). On the opposite, high delay solutions are located close together as well, but source-relay links are weaker, reducing end to end transmission probability.

Conclusion

The work presented in this paper has derived a flexible framework for evaluating the performance of a wireless ad hoc network with respect to several performance criteria. It has been designed to account for the broadcast nature of wireless communications and for an accurate interference characterization for the network. Our framework allows for the determination of both Pareto bounds and corresponding Pareto solutions. Network model and bounds for 1-relay and 2-relay networks have been assessed though simulations. The derived MO problem relies on an extended search space whose complexity mainly grows with N . Further works will consider proper choice for N and improved optimization approaches that better leverage the problem structure. However, compared to similar studies in mesh networks [5], the complexity of the proposed framework is scalable with respect to the number of flows in the network thanks to our definition of forwarding probabilities and our interference model.

The tightness (or pessimism) of the upper MO bound presented in this paper will be assessed in future works. Therefore, we are currently investigating combined fountain and network coding strategies to provide feasible transmission strategies which constitute lower MO bounds on performance.

use section* for acknowledgement

Figure 1 :

 1 Figure 1: Network model

Figure 2 :

 2 Figure 2: 1-Relay network

RD Inria 2 -Figure 3 : 2 -

 232 Figure 3: 2-Relay network

  s) the incoming flow rate received for each slot by each destination D d ∈ D.

  For l = N.|T | and m = |D||T |, Q is a non-zero l-by-l matrix, D is a non-zero l-by-m matrix, I is a m-by-m identity matrix and 0 a m-by-l zero matrix. Q is referred to as the relaying matrix and D as the arrival matrix.

  by-|T | diagonal matrix whose diagonal elements D u iD d give the probabilities for a packet transmitted by a node i in time slot u to arrive at destination D d . Consequently, D u iD d = p u iD d . The source matrix M S calculates the initial flow vector F (1) knowing the emission rate vector of any source in O. An initial flow vector of source S o ∈ O is obtained from M S with: F (1) = τ So • M S where M S is a |T |-by-(N + |D|).|T | transition matrix defined as M S = [Q S D S ] with Q S and D S the relaying and arrival matrix for the packets sent by the source, respectively. We have Q S = Q So1 • • • Q SoN and D S = D SoD1 • • • D SoD |D| where Q Soi follows the pattern given by (8) and D SoD d is a |T |-by-|T | diagonal matrix whose diagonal elements are D SoD d = p u SoD d . The multiplication of τ So with Q S computes the outgoing flows at the relays for source S o , i.e. it computes the probability for any relay to forward a packet directly received from S o . The multiplication of τ So with D S computes the flows coming into the destinations for source S o , i.e. it computes the probability for any destination to receive a packet directly from S o .

  1) • D I where I is a |D||T |-by-|D||T | identity matrix. The steady state flows transmitted by the relays -→ F R (∞) are directly obtained from the fundamental matrix. The steady state flows coming into the destinations ← -F D (∞) are obtained in two steps: propagating the steady state relay flows to the destination using D and adding to it the direct transmissions from the source to the destinations ← -F D (1) of time epoch 1. An element of ← --F u D d (∞) gives the steady state sum of all packets being received at destination D d in time slot u.

Inria

  of Pareto Front for 2-relay Analysis Result of Pareto Front for 2-relay

Figure 4 : 1 -

 41 Figure 4: 1-(top) and 2-relay (bottom) bounds.

Table 1 :

 1

	Main notations in this section Flow rate vector at time epoch s i (s) Outgoing flow rate vector after s time epochs F (s) -→ F u ← -F u i (s) Incoming flow rate vector after s time epochs
	M	Transition matrix
	Q	Relaying matrix
	M F	Fundamental matrix
	M S , D Source and arrival matrix
	4 Steady state evaluation

  1 .

	2-relay network The 2-relay network of Fig. 3 is considered again. The
	source is transmitting in time slot 1, relay A and B are transmitting in time slot
	2 and 3, respectively. The only non-null forwarding probabilities are x 23 AB , x 32 BA ,
	x 12 SA et x 13 SB . Initial flow vector F (1) is composed of
	Inria

Table 2 :

 2 RMSE for f * D and f *

	E

RR n°7799

Acknowledgments

This work was supported in part by the Marie Curie OIF Action of the European Community's Sixth Framework Program (DistMO4WNet project). This article only reflects the authors' views and the European Community are liable for any use that may be made of the information contained herein.


Carolina Bernal

Frédéric Christophoul
email: frederic.christophoul@get.obs-mip.fr

José Darrozes
email: darrozes@get.obs-mip.fr

Alain Laraque
email: alain.laraque@ird.frj.darrozes

Luc Bourrel
email: bourrel@get.obs-mip.fr

Jean-Claude Soula
email: soula@get.obs-mip.fr

Jean-Loup L Guyot
email: jean-loup.guyot@ird.fr

Patrice Baby
email: baby@ird.fr

Jean-Claude Soula

J-C Soula

Crevassing and

Keywords: Crevasse-splay, avulsion, anastomosis, Amazonian basin, floodplain, Pastaza

-INTRODUCTION

Avulsion is the process whereby a channel belt shifts relatively abruptly from one location in the floodplain surface to another in favour of a steeper gradient [START_REF] Slingerland | River avulsion and their deposits[END_REF][START_REF] Mohrig | Interpreting avulsion process from ancient alluvial sequences: Guadalope-Matarranya (Northern Spain) and Wasatch Formation (western Colorado)[END_REF][START_REF] Bridge | Rivers and Floodplains: Forms, Processes and Sedimentary Record[END_REF][START_REF] Stouthamer | Avulsion: The relative roles of autogenic and allogenic processes[END_REF]. Avulsions may result in a new channel belt [START_REF] Makaske | Anastomosing rivers: a review of their classification, origin and sementary products[END_REF].

Avulsions can take place in a variety of alluvial environments, including alluvial fans, braided river floodplains, meandering river floodplains, anastomosed rivers and delta plains [START_REF] Törnqvist | Middle and Late Holocene avulsion history of the River Rhine (Rhine-Meuse delta, Netherlands)[END_REF][START_REF] Jones | Channel avulsions and related processes, and large -scale sedimentation patterns, Rio Grande, San Luis Valley, Colorado[END_REF][START_REF] Morozova | Holocene avulsion style and sedimentation patterns of the Saskatchewan River, Cumberland Marshes, Canada[END_REF][START_REF] Makaske | Anastomosing rivers: a review of their classification, origin and sementary products[END_REF]Jain & Sinha, 2005;[START_REF] Leier | Mountains, monsoons and megafans[END_REF][START_REF] Assine | River avulsions on the Taquari megafan, Pantanal wetland, Brazil[END_REF][START_REF] Sarma | Fluvial process and morphology of the Brahmaputra River in Assam, India[END_REF][START_REF] Sinha | Geomorphic characterization and diversity of the rivers of the Gangetic plains[END_REF]. Avulsions are primarily features of aggrading floodplains [START_REF] Slingerland | River avulsion and their deposits[END_REF] where they constitute the main cause of anastomosing channel initiation [START_REF] Makaske | Anastomosing rivers: a review of their classification, origin and sementary products[END_REF].

River avulsion setup is a long term process where the river aggrades over tens to thousands of years and becomes poised for avulsion. Several mechanisms have been proposed as being responsible for important aggradation rates including increased sediment load relative to capacity, change in peak water discharge, or decrease in channel gradient owing to sinuosity increase, delta extension, base level perturbation, downstream tectonic uplift [START_REF] Jones | Causes of avulsion: an overview[END_REF][START_REF] Stouthamer | Factors Controlling the Holocene Avulsion History of The Rhine-Meuse Delta (The Netherlands)[END_REF][START_REF] Slingerland | River avulsion and their deposits[END_REF]. The trigger is a short-term event which immediately provokes the avulsion. Triggers are most frequently floods, but could also be processes as varied as abrupt tectonic movements, ice jams, log dams; vegetative blockages, debris dam, beaver dams, hippopotami trails, or other bank failures and downstream migration of bars that temporarily block the throat of a branch [START_REF] Stanistreet | Low sinuosity meandering river bedload rivers of the Okavongo Fan: channel confinment by vegetated levées without fine sediment[END_REF][START_REF] Jones | Causes of avulsion: an overview[END_REF][START_REF] Mohrig | Interpreting avulsion process from ancient alluvial sequences: Guadalope-Matarranya (Northern Spain) and Wasatch Formation (western Colorado)[END_REF][START_REF] Slingerland | River avulsion and their deposits[END_REF].

It has been suggested that two end-members of the avulsion series could be (1) incisional avulsion characterized by an early incision and then infilling by migrating bar and (2) aggradational avulsion beginning with aggradation followed by integration of streams into a single downcutting channe [START_REF] Smith | Anatomy of an avulsion[END_REF]. [START_REF] Mohrig | Interpreting avulsion process from ancient alluvial sequences: Guadalope-Matarranya (Northern Spain) and Wasatch Formation (western Colorado)[END_REF] proposed that the type of avulsion is strongly influenced by whether or not the adjacent flood plain is well drained.

Several numerical models of avulsion have been proposed in order to understand mechanics and sedimentary architecture of alluvial plains, deltas, and fans (Slingerland and Smith, 1998;[START_REF] Jerolmack | Complexity in a cellular model of river avulsion[END_REF]. However, these models do not include some aspects of the fluvial environment (e.g., vegetated levees and floodplains) that may play a key role in avulsion dynamics [START_REF] Mohrig | Interpreting avulsion process from ancient alluvial sequences: Guadalope-Matarranya (Northern Spain) and Wasatch Formation (western Colorado)[END_REF]. Indeed avulsions although infrequent [START_REF] Smith | Anatomy of an avulsion[END_REF][START_REF] Morozova | Holocene avulsion style and sedimentation patterns of the Saskatchewan River, Cumberland Marshes, Canada[END_REF] can best be studied by direct observation [START_REF] Törnqvist | Middle and Late Holocene avulsion history of the River Rhine (Rhine-Meuse delta, Netherlands)[END_REF][START_REF] Mohrig | Interpreting avulsion process from ancient alluvial sequences: Guadalope-Matarranya (Northern Spain) and Wasatch Formation (western Colorado)[END_REF].

The present study focuses on the lower Río Pastaza (Ecuadorian and northern Peruvian Amazonia) which drains one of the largest tropical alluvial megafan: the Pastaza Megafan [START_REF] Räsänen | Evolution of the Western Amazon Lowland Relief : impact of Andean foreland dynamics[END_REF][START_REF] Bernal | Late glacial and holocene avulsions of the Rio Pastaza Megafan (Ecuador-Peru): Frequency and controlling factors[END_REF]. Recently, substantial attention has been paid to the construction of the Pastaza Megafan (Christophoul et al., 2002;Bes de Berc et al., 2005;[START_REF] Burgos | Genese et progradation d'un cone alluvial au front d'une chaine active: example des Andes Equatoriennes au Neogene[END_REF][START_REF] Bernal | Late glacial and holocene avulsions of the Rio Pastaza Megafan (Ecuador-Peru): Frequency and controlling factors[END_REF]. In the present study, an attempt is made to depict recent avulsions in the Lower Rio Pastaza, near its confluence with the Rio Marañon and to discuss their nature and their role in the construction of floodplain. This study is based on the analysis of the evolution of the fluvial morphology using satellital images and supported by the hydro-pluviometric data yet recorded in the region.

-GEOLOGIC AND GEOMORPHOLOGIC CONTEXT

The Río Pastaza drains the eastern Ecuadorian Andes and the upper Amazonian Basin (Fig. 1). The Ecuadorian Andes are characterized by recently active volcanoes, two of them, the Tungurahua (5,020 m a.s.l.) and the Sangay (5,230 m a.s.l.) having had eruptions all through the 20 th century. Fluvial aggradation started in the Amazonian foreland around 22 Ma with sediments derived from the west (Christophoul et al., 2002[START_REF] Roddaz | Forebulge dynamics and environmental control in Western Amazonia: The case study of the Arch of Iquitos (Peru)[END_REF][START_REF] Roddaz | Cenozoic sedimentary evolution of the amazonian foreland basin system[END_REF][START_REF] Burgos | Genese et progradation d'un cone alluvial au front d'une chaine active: example des Andes Equatoriennes au Neogene[END_REF]. This sedimentation corresponds to the deposits of the Neogene Napo-Pastaza megafan sensu lato. The megafan apex migrated forward as a result of thrust fault propagation (Bès de [START_REF] De Berc | Geomorphic evidence of active deformation and uplift in a modern continental wedge-top-foredeep transition: Example of the eastern Ecuadorian Andes[END_REF] but no tectonic influence has been detected in the study area except flexural subsidence origined as a result of thrust stacking in the cordilleras [START_REF] Roddaz | Forebulge dynamics and environmental control in Western Amazonia: The case study of the Arch of Iquitos (Peru)[END_REF]. The most recently active part of the Pastaza megafan [START_REF] Räsänen | Evolution of the Western Amazon Lowland Relief : impact of Andean foreland dynamics[END_REF][START_REF] Räsänen | Recent and ancient fluvial deposition systems in the Amazonian foreland basin[END_REF][START_REF] Zapata | Evolución del Abanico del Pastaza y su relación con los depósitos Ucamara[END_REF][START_REF] Bernal | Late glacial and holocene avulsions of the Rio Pastaza Megafan (Ecuador-Peru): Frequency and controlling factors[END_REF] is located in the southeast of the Neogene megafan, in the subsiding area between the Pastaza -Marañon -Ucayali confluence and the Iquitos Arch. The Rio Pastaza is 710 Km long and its drainage basin covers 39,000 Km 2 . The modern course of the Rio Pastaza in the Amazonian plain is the result of numerous avulsions (more than 100) which occurred since the Last Glacial Maximum [START_REF] Bernal | Late glacial and holocene avulsions of the Rio Pastaza Megafan (Ecuador-Peru): Frequency and controlling factors[END_REF]. The last major avulsion of the Rio Pastaza in the Amazonian plain is older than AD 1691 [START_REF] Bernal | Late glacial and holocene avulsions of the Rio Pastaza Megafan (Ecuador-Peru): Frequency and controlling factors[END_REF].

-THE RIO PASTAZA BASIN

In the Amazonian foreland, the mean slope of the Rio Pastaza progressively decreases from 2.0 m / km at its exits from the cordillera to 0,2 m / km at its confluence with the Rio Marañon (Fig. 2). In the north, the Rio Pastaza forms a long (125 km) and smooth southward curve and the planform morphology changes from braided to low sinuosity single channel. (Fig. 2, reaches C and D).

In the south (Fig. 2, reach E), the Rio Pastaza preserves a main channel, a few metres higher than the flood plain, and several diverging single channels bifurcate from this main channel on both banks to re-join it downstream thus forming a typical anastomosed system.

-HYDROMETEOROLOGICAL DATA

A previous study of the pluviometry and discharge in Ecuadorian streams [START_REF] Laraque | Sediment budget of the Napo river (Amazon basin, Ecuador and Peru)[END_REF] has shown that rainfall and discharge regimes are highly variable in space and show a strong orographical control. The pluviometric and hydrologic regimes of the Río Pastaza in various parts of its drainage basin have been established using data issued from three pluviometric and two hydrosedimentological stations (location on Fig. 2) over the period . The data issued from the station situated in the range (Cusubamba and Baños pluviometric stations, Baños hydrologic station) are highly dependant upon their situation. At the proximity of the study area, it is practically impossible to get any data because of isolation and insecurity [START_REF] Laraque | Sediment budget of the Napo river (Amazon basin, Ecuador and Peru)[END_REF]. In contrast, the stations of the middle Pastaza basin also situated in the Amazonian wetland can be considered to be representative of the lower Pastaza basin.

At the Teniente Lopez meteorological station (230 m a.s.l.), the average annual rainfall is 2,579 mm and is distributed roughly evenly throughout the year (7-10% each month, Fig. 3). The minimum and maximum values of average monthly rainfall are 188.6 mm (March) and 252.8 mm (April), respectively [START_REF] Espinoza | Spatiotemporal rainfall variability in the Amazon Basin Countries (Brazil, Peru, Bolivia, Colombia and Ecuador)[END_REF]. At the Andoas hydrological station (Fig. 2) 150 km upstream of the study area sediment discharge has been recorded for 5 years (1997)(1998)(1999)(2000)(2001)(2002), with an average annual discharge of 910 m 3 s -1 . The highest average monthly discharge (Fig. 3) occurs in May (1,222 m 3 s -1 ) and the lowest in September (674 m 3 s -1 ).

Remote sensing analyses of the hydrology of the Amazonian wetland show that the lower Pastaza is like the rest of Amazonia affected by annual floods invading the whole floodplain. At high water stage, the floodplain entire of distal part of the Pastaza Megafan and as consequence, the lowermost Pastaza is covered by water [START_REF] Mellack | Remote sensing of wetlands on a global scale[END_REF][START_REF] Martinez | Mapping of flood dynamics and spatial distribution of vegetation in the Amazon floodplain using multitemporal SAR data[END_REF]. In contrast with the range and the immediate piedmont, we have no data showing the impact of ENSO events.

-DATA AND METHODS

-Satellite images

This study was based on the integration of available hydro-sedimentological and geological information with new morphologic and topographic data obtained from satellite and/or aerial images and a digital elevation model.

Spatial variations of channel pattern as well as avulsion occurrences were traced over the past 31 years using ASTER and Landsat time series satellite imagery (MSS, TM, ETM+, 1977-2002), and satellites CBERS (China-Brazil Earth Resources Satellite, 2008) for Ecuador-Peru. We selected the less cloudy images for approximately the same hydrological period (except the 2002 satellite image).

When necessary, the selected images were orthorectified using the software ENVI 4.2; and projected to UTM Zone 18 S. The channel's vectors were digitized using the Geographic Information System (Arc View 9.1). However one of the main problems was the discrimination between open water and land. In order to solve this we used the wavelengths in the near and middle-infrared ranges (740 -2500 nm; [START_REF] Jensen | Remote Sensing of the Environment: An Earth Resource Respective[END_REF]. However, the Andean piedmont remains very cloudy at any time and optical/IR image acquisitions remain uncertain, and the temporal coverage is weak. This supports the systematic use of SRTM images in spite their lower resolution.

-Topographic data

Topographic data used in this study are those of the SRTM-3 V4 DEM [START_REF] Jarvis | Hole-filled seamless SRTM data V4[END_REF]. Interest of SRTM-3 data for geomorphological studies even in low relief areas have been shown by [START_REF] Coltelli | SIR-C/X-SAR multifrequency multipass interferometry: A new tool for geological interpretation[END_REF], [START_REF] Farr | The Shuttle Radar Topography Mission[END_REF] and Valeriano et al., (2007). These data are based on radar interferometry (INSAR) in C band. Pixel size at the equator is 90m and slightly increases with the latitude. According to [START_REF] Rodriguez | An assessment of the SRTM topographic products[END_REF], error in height, in the studied part of the Amazonian lowland is ~ 5m.

The SRTM DEM for our study area was reprocessed with ENVI 4.3 using elevation issued from topographic maps of the Ecuadorian and Peruvian geographic surveys. DEM SRTM V4 data are shown on Figure 4. In order to check the coherence of elevation data provided by the SRTMV4 DEM, we extracted a river long profile along the left bank of the main channel of the lowermost Pastaza (Fig. 5). Indeed, open water bodies provoke a loss of coherence in radar response and their elevation in the DEM results from interpolation with surrounding data. In the case of an alluvial river such as the lower Pastaza it should be reasonably assumed that river bank slope and river slope are similar, even though this is not necessarily the case in incising rivers. In fact, the slope computed from the elevation data is 0.25 m.km -1 (Fig. 4) and that extracted from topographic maps is 0.2 m.km -1 (Fig. 2).

-IMAGE AND TOPOGRAPHY ANALYSES

Comparing satellite images allowed us to define three areas of interest along the lower Pastaza (Fig. 2 reach E). The first example is located in the right bank of the Rio Pastaza ( 1 on Fig. 5) and is concerned with opening of a channel joining the Rio Pastaza and one of its tributaries, the Rio Huangana. The second example is located downstream, to the left of the Rio Pastaza (2 on Fig. 5), and illustrates the evolution of a crevasse-splay into an anastomosed channel. The third example is an ancient crevasse splay more evolved than the preceding having given rise to a larger-scale anabranch downstream of the formers (3 on Fig. 5).

Despite the scarcity of cloud-free images 10 scenes showing a marked change in the fluvial morphology were found for the period 1977-2008 in the first example (Rio Pastaza-Rio Huangana junction) and 8 for the second. No optical images showing morphological changes have been found for the third example but this lack was supplied by a topographic analysis allowed by the SRTM DEM.

6.1-The Rio Huangana crevasse-splay

The initial image (1981/07/21) (Fig. 6) shows the Río Pastaza and the Río Huangana flowing almost parallel to each other. The Río Huangana is separated from the Pastaza by a large natural levee but will rejoin the Pastaza downstream. 1987/07/18 and 1989/07/23 images (Fig. 6) show a progressive diminution in the water flow of the Rio Huangana (Fig. 6). The July 1990 image (Fig. 6) shows a small diverging stream "A" issued from the outer bend of the Rio Pastaza and breaching the levee to join the Rio Huangana downstream.

The 1991, 1992 & 1994 images (Fig. 6) show a slow but constant widening of the stream A. As a consequence, there is a considerable increase in water volume inducing the development of a crevasse-splay by the "A" channel (2001/06/30, Fig. 6) at the confluence with the Huangana, associated with the creation of a marsh between the Pastaza's western levee and the floodplain.. This marsh may have been a result of the formation of the crevasse splay damming the upper reach of the Huangana already situated at a lesser elevation than the Pastaza main channel). The 2008/02/01 image shows the Huangana River beheaded (Fig. 6). The distributive channels forming the crevasse-splay evolved into an anastomosing network rejoining the Rio Pastaza along the downstream course of the Rio Huangana.

-Splay 2, left bank of the Rio Pastaza

Figure 5.2 shows a 26.16 km 2 large crevasse-splay made up of 3 lobes formed from 3 alluvial ridges. The 1981/07/21 image (Fig. 7) shows a well-developed crevasse-splay with a distributive network prograding into the floodplain. The July 1987, 1899 and 1990 images (Fig. 7), show these undefined flooded areas have disappeared and the trunk channel bifurcating downstream, leading to active channels organized in a classical distributary pattern. Only in 1990 (Fig. 7) the trunk channel joins the Pastaza. Between 1991 and 2008 the network becomes simpler with the progressive consolidation of a single trunk channel, which joins the Pastaza River downstream (Fig. 7).

The topography of this crevasse-splay is shown in figure 4. 4 transverse and 4 river/alluvial ridge long profiles were extracted from the DEM (Fig. 8). Profile 1 (Fig. 8) is close to the apex of the crevasse-splay and reveals a 1.7 km wide alluvial ridge culminating at 158 m a.s.l. above a low relief area with a mean elevation between 142 and 147m. Profile 2 (Fig. 8) is located in the area where the initial alluvial starts dividing into 2 and then 3 second order alluvial ridges. The ridge as a whole is now 2.9 km wide and culminates at 155 m a.s.l.. Profile 3 (Fig. 8) exhibits the 3 second order ridges. Each of these ridges is around 1km wide and culminate at 155 m a.s.l. Profile 4 crosses the two south-western lobes of the crevasse-splay which are wider (2 and 3.2 km) and less elevated (147 and 149 m a.s.l.) than the upstream ridges whereas the floodplain is 136 m a.s.l. high. A V-shaped depression is observed between the eastern second order ridges. This depression is interpreted as the non-vegetated active channel (lighter green colour in the Landsat image, Fig. 5.2). A smaller but similar depression appears in the western second order ridge (profile 4,Fig. 8,left). This would indicate that the channels widen downstream within the alluvial ridges.

Profiles Lp1 to Lp4 are long profiles of the alluvial ridges (Fig. 8). The upper enveloping line of the western more profile Lp1 shows a gentle downstream slope. This shape is globally concave up with a slope increasing downstream. The mean slope is 1.64.10 -3 that is an order of magnitude higher than the slope of the Rio Pastaza (0.2.10 -3 , see fig. 5A). Profile Lp2 is that of the central lobe of the crevasse-splay (Fig. 5). Its overall shape is similar to that of Lp1 and its mean slope is 1.22.10 -3 . Profile Lp3 is that of the eastern lobe of the crevasse-splay. It is shorter than the other alluvial ridges of the splay (2.3 km) and has a steeper slope (2.07.10 -3 ).

-Splay 3, southern left bank of the Rio Pastaza

This crevasse-splay (fig. 5.3) and its structure are difficult to identify on satellite images because of its low contrast, which suggests that vegetation is quite homogeneous and higher (darker) all over the splay. In contrast, the SRTM DEM reveals very clearly the structure (Fig. 4). 4 transverse topographic profiles (Profiles 5 to 8) and ridge long profiles (Lp4 to Lp8) are shown in Figure 9. Profiles 5 to 8 depict a proximal-distal gradient within the crevasse-splay. Profile 5 exhibits the same bell shape than profile 1 of Splay 2 and can similarly be interpreted as an alluvial ridge. The top of this alluvial ridge culminates at 157 m whereas the floodplain has an elevation of 136 m. Profile 6, located downstream in the crevasse-splay, enables the identification of two alluvial ridges. The first one is 5 km wide and culminates at 155 m. The second is 1.7 km wide. Profile 7 is located still downstream and crosses only one alluvial ridge 2 km wide and 16 m high). The same alluvial ridge 1.2 km wide and 12 m high is visible in Profile 8 where another ridge is observed. The alluvial ridge long profiles (Lp4 to Lp8 on fig. 9) reveal the morphological characteristics of these splays. Profiles Lp4 and Lp5 (Fig. 9) follow short transverse ridges (~3.5Km) formed in the northern part of the splay (Fig. 5). These ridges are directly connected to the crevasse channel through the right bank levees of the Rio Pastaza main channel. They exhibit regularly decreasing slopes (4.85.10 -3 and 6.28.10 -3 , respectively) toward the floodplain. This slope is an order of magnitude higher than the Pastaza river gradient (0.25.10 -3 , Fig. 5). Lp6 (Fig. 9) shows the longest (longitudinal) alluvial ridge of the splay. The first 12 km of the profile exhibit a low but constant decrease in elevation with a 0.66.10 -3 slope followed by an abrupt decrease in elevation in the last kilometre. This decrease in elevation is interpreted as a vegetation change from the alluvial ridge to the floodplain. Lp7 and Lp8 are along longitudinal ridges/channels bifurcating from the Lp6 ridge/channel in the southern lobes of the splay. Lp7 is 3.8 km long. Its slope is 3.27.10 - 3 . Lp8 is 5.5km long and has a 1.24.10 -3 slope.

Satellite images reveal no evolution of the distributive drainage network of the crevasse-splay through time. This apparent stability of the drainage network along with the homogeneous vegetation suggest that this crevasse-splay is older than the other crevasse-splays described above and represents the terminal stage of evolution of a crevassesplay into an anastomosed channel. The older age of this crevasse splay may explain why Lp6 ridge/channel (fig. 9) has a lower slope than the other ridge/channels. It suggests that the anastomosed channel tries to reach an equilibrium profile with a regular slope by smoothing the knickpoint initially created at the jonction of the crevasse-splay distributary and the floodplain drain it annexes.

-DISCUSSION

8-1 -History of avulsions

All the 3 avulsion sets analysed involve crevasse splays showing the same evolutional trend: (1) initial crevassing; (2) crevasse splay propagation in the floodplain; (3) capture of one of the splay channels by a pre-existing channel of the floodplain (or annexation) with headward erosion progressively reducing the initial channel slope. Because headward incision diminishes the elevation of the captured splay channel, the flow of the other distributaries of the splay is captured into this channel (Fig. 10) (see also [START_REF] Smith | Anatomy of an avulsion[END_REF][START_REF] Slingerland | River avulsion and their deposits[END_REF]. In our study area, this process will then arrive to the formation of a diverging -rejoining channel being part of a larger-scale anastomosing system.

In detail, the long profiles of the alluvial ridges indicate that an alluvial ridge evolved into an anastomosing channel after having acquired a lower slope than the others (Lp6 on Fig. 9 and in a lesser extend in Lp2 on Fig. 8). [START_REF] Slingerland | River avulsion and their deposits[END_REF] propose that gradient decrease by headward erosion be associated with avulsions by incision and did not mention (but did not reject) this process for avulsions by progradation and annexation. In our cases, the difference in elevation between the floodplain and the ridge containing the avulsed channel will produce profile readjustment in any case. In fact, the systems studied experienced incision, annexation and progradation at various stages in their evolutional history. Splay 1 (Rio Huangana) started as an incising channel, then annexing a pre-existing river, and finally forming a prograding splay (beheading the host river) as a result of increasing flow and sedimentary load. Splay 2 initiates as a progradational splay with one of the distributive channels being then captured by (or annexing) a pre-existing small stream of the floodplain, implying incisional profile readjustment. Splay 3 likely initiated as a progradational splay, with a channel captured by (or annexing) a stream of the floodplain then affected by a strong incision.

-Origin of avulsions

Set up and cause. The lower Pastaza flows in an area affected by flexural subsidence where aggradation rate is important. As for other highly aggrading alluvial rivers, the lower Pastaza main channel appears to be surelevated with respect to the floodplain, inducing a high ratio between channel-to-floodplain slope and channel slope (see [START_REF] Jones | Causes of avulsion: an overview[END_REF]. This channel is bounded by poorly vegetated levees. In the 3 cases we have observed that breaching of the levees occurred in the outer arc of meanders where the shear stress at the top of the levvees during flodds is maximum. This location seems to be at the origin of the rather regular anastomoses of the lower Pastaza on both sides of the main channel.

Trigger. In the studied area the seismic activity is reduced, log dams, vegetative blockages, animal trail are unlikely because of the size and high stream power of the river. Seismic fault surface rupture are also unlikely to occure in an area located far eastward from the subandean active faults. Floods are therefore the most likely triggers of avulsions. In our study area there are no evidences of exceptional floods related or not to ENSO events. In contrast, it has been shown that annual floods occur in this area as in the rest of the Amazonian wetland. It is therefore likely that annual floods triggered the avulsions. As reported by [START_REF] Jones | Causes of avulsion: an overview[END_REF], most of the observed avulsions are triggered by 'ordinary' floods provided that channel instability was close to the avulsion threshold. Meander outer arcs will constitute weakness zones in which crevassing and then avulsion will most easily occur. [START_REF] Aalto | Episodic sediment accumulation on Amazonian floodplains influenced by El Niño/Southern Oscillation[END_REF] pointed out the role of ENSO in the Amazonian Basin bringing great amount of water and sediment in amazonian rivers of Bolivia. Data presented in this sudy prevent us to correlate vrevassing and ENSO Event. Indeed, discharge data in Andoas hydrologic station do not cover the last significative ENSO Events. Moreover, satellite images cannot help us to show initial crevassing is related to ENSO events. The only case in which such event could have influenced the evolution of a crevasse splay would be the Rio Huangana. The development of the splay in the floodplain occurred between 1994 and 2001 (Fig. 6). During this timerange, a sevre ENSO occurred in in the winter 1997-98.

-CONCLUSION

In spite of minor differences, the three examples studied in the present paper show the same evolutional history: (1) initial crevassing; (2) crevasse splay propagation in the floodplain; (3) capture of one of the splay channels by a pre-existing channel of the floodplain with headward erosion progressively reducing the initial channel slope; (4) change of this channel into a diverging -rejoining channel being part of a larger-scale anastomosing system. The main difference is concerned with the capturing (or annexed) floodplain stream. In example 1 (Rio Huangana) this stream is a tributary flowing parallel to the main channel for a rather long distance before being annexed and beheaded. The splay is rather distant from the main channel, relatively small, and immediately postdates the annexation. In examples 2 and 3, the splay forms where crevassing occur, is relatively large, and develops a distributary network. Capture or annexation selects one of these drains to further become an anabranch of the future larger-scale anastomosing system. Crevassing and consecutive avulsions take advantage of the presence of meanders and this location seems to be at the origin of the rather regular anastomoses of the lower Pastaza.

Although the precipitation regime and discharges curves are rather regular, levees breaching and crevasse splay formation are likely to be results of annual floods which are present in this area as in the rest of Amazonia.

Figure 1 -

 1 Figure 1 -Mosaic of Digital Numeric Model of Napo-Pastaza Megafan. The watershed is represented in black lines and the modern river network, in white. White squares represent active volcanoes (Reventador Volcano, RV; Tungurahua Volcano, TV; Sangay Volcano, SV) and blue spots represents active volcanoes. Black dashed line: boundary of the rio pastaza river basin. Red dashed line is the boundary of the Napo-Pastaza megafan (Bes de Berc et al., 2005). Yellow dashed line represents the post-LGM Pastaza Megafan[START_REF] Bernal | Late glacial and holocene avulsions of the Rio Pastaza Megafan (Ecuador-Peru): Frequency and controlling factors[END_REF].

Figure 2 -

 2 Figure 2 -The Rio Pastaza: drainage basin, hypsometric profile of the Rio Pastaza and associated channel patterns (reaches A to E). White squares:active volcanoes with TV: Tungurahua and SV: Sangay. Black crosses:(Cusubamba St., Cu; Baños St., Ba; Andoas St., An; Teniente López St., TL).

Figure 4 -

 4 Figure 4 -SRTM V4 DEM of the lower Rio Pastaza and topographic profiles. This DEM shows the topography of the Modern Rio Pastaza alluvial ridge and its associated floodplain. On the left bank of the Rio Pastaza, two crevasse-splays can be observed. The correspond respectively to the crevasse-splay 2 and 3 on Figure 3. The crevasse-splays are bounded by a light blue dashed lines. Yellow solid lines 1 to 8 and black solid lines Lp1 to Lp6 corresponds to the topographic profiles exhibited on figures 8 and 9. A and B: topographic profiles with A: topographic profile across the Rio Pastaza alluvial ridge, and B: topographic profile of the Rio Pastaza alluvial ridge extracted along the Rio Pastaza right bank.

Figure 5 -

 5 Figure 5 -Location and satellite images of the 3 crevasses splays and anastomosing channels detailed in this study (Landsat TM, VNIR band). Black dashed line: boundary of the Rio Napo drainage basin, green solid line: active channels. 1) Rio Huangana crevasse-splay, 2) Splay 2, left bank of the Rio Pastaza, 3) Splay 3, Southern left bank of the rio Pastza. Yellow dashed line: boundary of crevasse-splays. CS: crevasse-splay, CCh: Crevasse Channel, AvC: Anastomosing channel.

Figure 6 -

 6 Figure 6 -Time series of satellite images (Band 5 -Landsat & Band 4 -CBER) recording the avulsion and crevasse-splay evolution of the Rio Huangana crevasse-splay on the right bank of the Rio Pastaza. On the left: map of the lowermost Rio Pastaza, exhibiting the Rio Huangana, its point of divergence with the Rio Pastaza (red box) illustrated by the time series on the right. Remark the downstream junction of the Rio Huangana with the Rio Pastaza through the Rimachi floodplain lake.

Figure 7 -

 7 Figure 7 -Time series of satellite images (Band 5 -Landsat & Band 4 -CBER) recording the avulsion evolution and crevasse-splay evolution of the Splay 2 on the northern left bank of the Rio Pastaza.

Figure 8 -

 8 Figure 8 -Topographic profiles of the Splay 2 on the northern left bank crevasse-splay. Location of profiles on figure 5. Profiles 1 to 4: topography across the crevasse-splay. Profiles Lp1 to Lp3: river long profiles of the alluvial ridge of the crevasse-splay.

Figure 9 -

 9 Figure 9 -Topographic profiles of the Splay 3, on the southern left bank crevasse-splay. Location of profiles on figure 5. Profiles 5 to 8: topography across the crevasse-splay. Profiles Lp4 to Lp8: river long profiles of the alluvial ridge of the crevasse-splay.

Figure 10 -

 10 Figure 10 -Sketch map showing the evolution of the lower Pastaza crevasse-splay into an anastomosing alluvial ridge. AR: Alluvial ridge (including channel and natural levees), CS: crevasse-splay, DN: Distributary network, FD: Floodplain drainage. A: initial crevassing and progradation of the crevasse-splay in the floddplain. B: crevasse-splay channel capture by a pre-existing floodplain drain and abandonment of the crevase splay channel after the capture. C: final abandonment of the crevasse-splay channels end development of a new alluvial ridge joining the former crevasse-splay drain upstream and the former floodplain drain downstream.

Acknowledgments

This investigation was supported by a doctoral Albαn fellowship of the European Union (n° E05D057404EC), IRD and Université Paul Sabatier grants. We are grateful to the HYBAM/Ecuador project and its partnership institutions, INAMHI and IRD; and to the 'Instituto Geofisico' of the 'Escuela Politecnica Nacional' of Quito -Ecuador. The authors thank Centre for Tropical Agriculture (CIAT) for providing topographic data through http :// srtm . csi . cgiar . org


Jan Bartovský

Eva Dokládalová

Michel Bilodeau

Petr Dokládal

Texture Analysis with Arbitrarily Oriented Morphological Opening and Closing

Keywords: Mathematical Morphology, Granulometry, parallel implementation

This paper presents a fast, streaming algorithm for 1-D morphological opening on 2-D support. The algorithm is further extended to compute the complete size distribution during a single image run. The Structuring Element (SE) can be oriented under arbitrary angle that allows us to perform different orientation-involved image analysis, such as local angle extraction, directional granulometries, etc.

The algorithm processes an image in constant time irrespective of the SE orientation and size, with a minimal latency and very low memory requirements. Regardless the SE orientation, it reads and writes data strictly sequentially in the horizontal scan order. Aforementioned properties allow an efficient implementation in embedded hardware platforms that opens a new opportunity of a parallel computation, and consequently, a significant speed-up.

Introduction

During several last decades, the problems of object recognition, texture analysis, orientation, and size distribution have become very important in image analysis and computer vision field. Tthe Mathematical Morphology (MM) has proved to be a very powerful tool in solving such problems [START_REF] Matheron | Random sets and integral geometry[END_REF] [START_REF] Maragos | Pattern spectrum and multiscale shape representation[END_REF]. On the other hand, the computation complexity of profound analysis over many shapes and rotations is overwhelming as MM often resides on iterative use of two basic operations dilation and erosion.

The paper is organized as follows: Section 2 recalls the basic notion of morphological operators. Section 3 outlines the main principle of the proposed algorithm. Section 4 presents performance results of computer benchmarks and a comparison with other algorithms. Finally, Section 5 displays a few examples of opening and size spectrum applications.

State of the Art

One of the early algorithms of morphological erosion that complete in O(1) per pixel is van Herk [START_REF] Van Herk | A fast algorithm for local minimum and maximum filters on rectangular and octagonal kernels[END_REF]. Although it runs in constant time using only 3 comparisons per pixel, it requires two passes: direct and reverse, which degrades its real-time features especially if extended to 2-D. Later, Gil and Kimmel [START_REF] Gil | Efficient dilation, erosion, opening, and closing algorithms[END_REF] improved van Herk algorithm and reduced the number of comparison to only 1.5 per pixel. Soille et al. [START_REF] Soille | Recursive implementation of erosions and dilations along discrete lines at arbitrary angles[END_REF] first published the principle of arbitrary oriented SE. The SE rotation is handled by Bresenham algorithm, the computation is carried out by the van Herk algorithm.

Van Droogenbroeck and Buckley [START_REF] Van Droogenbroeck | Morphological erosions and openings: Fast algorithms based on anchors[END_REF] proposed a different approach based on notion of anchors, the points that are not affected by a transformation. This algorithm uses a histogram for calculation making any implementation dependent on the number of gray levels.

Urbach and Wilkinson [START_REF] Urbach | Efficient 2-d grayscale morphological transformations with arbitrary flat structuring elements[END_REF] proposed an algorithm for arbitrary shaped 2-D flat SE based on the computation of multiple horizontal line SEs for every pixel and storing them in a lookup table. The result is then computed by taking the maximum from the lookup table corresponding to the shape of the SE.

The issue of granulometry has been already targeted by Vincent [START_REF] Vincent | Granulometries and opening trees[END_REF] and Menotti-Gomes et al. [START_REF] Menotti-Gomes | 1D Component tree in linear time and space and its application to gray-level image multithresholding[END_REF]. Both computed the size distribution and opening trees in linear time with respect to the size of the image.

Novelty

We propose a new algorithm that computes both the opening and size spectrum at the same time with the following properties. The algorithm supports line SE oriented in arbitrary angle. It processes an image in constant time with respect to length of SE or its orientation. It accesses to the input and output data strictly sequentially in a horizontal scan order inferring a small computation latency. This derandomization of memory bandwidth along with very low memory requirements, which are much lower than the mere size of the image to process, make the algorithm especially interesting for realtime hardware systems. The single-thread C implementation of this algorithm brings a very low execution time and outperforms a few other existing efficient algorithms in the case of opening. In the case of size spectrum, the speed-up against the conventional approach, which requires many openings to be computed and then subtracted (4), is in orders of magnitude.

Basic Notions

Let δ B , ε B : Z 2 → R be a dilation and an erosion on gray-scale images, parameterized by a SE B assumed flat (i.e., B ⊂ Z 2 ) and translation-invariant, defined as

δ B (f ) = b∈B f b ; ε B (f ) = b∈ B f b . (1)
The hat denotes the transposition of the structuring element, equal to the set reflection B = {x | -x ∈ B}, and f b denotes the translation of the function f by some scalar b.

The SE B is equipped with an origin x ∈ B.

Let ϕ B , γ B : Z 2 → R be a closing and an opening on gray-scale images, parameterized by the aforementioned SE B, defined as

ϕ B (f ) = ε B [δ B (f )] ; γ B (f ) = δ B [ε B (f )] . (2)
Let S λB : R 2 → R be a size spectrum, parametrized by a SE B ⊂ R 2 and its size λ, defined as (originally proposed by Maragos [START_REF] Maragos | Pattern spectrum and multiscale shape representation[END_REF])

S λB (f ) = - d dλ γ λB f ; f : R 2 → R (3)
Since we are interested in images with a discrete and bounded support

D ⊂ Z 2 , D = [1..M ] × [1.
.N ], the discrete size spectrum S λB is transformed to as follows

S λB (f ) = D γ (λ-1)B f -γ λB f ; f : D → R, (4)
considering the size spectrum step dλ = 1.

Hereafter, we focus our description on the opening by a line SE B α λ . This SE has shape of a discrete line of length λ rotated by angle α from positive x-axis counterclockwise. The closing can be obtained in accordance to the duality property.

The oriented 1-D opening can also be used for detection of local orientation (orientation field) ζ λ (f ), see [START_REF] Soille | Recursive implementation of erosions and dilations along discrete lines at arbitrary angles[END_REF], image restoration, or oriented spectrum [START_REF] Van Droogenbroeck | Morphological erosions and openings: Fast algorithms based on anchors[END_REF]. The oriented size spectrum OS(α, λ) of an anisotropic texture is the size distribution expectancy of a 1-D signal obtained by intersection with a randomly drawn straight line. The expectancy is approximated by the frequency count in [START_REF] Gil | Efficient dilation, erosion, opening, and closing algorithms[END_REF].

ζ λ (f ) = arg max α∈(0,180) γ α λ (f ) (5)
[OS(α, λ)] (f ) = D γ B α λ-1 f -γ B α λ f f : D → R, (6)

Principle of the 1-D Opening Algorithm

We describe the main principles of our algorithm in this section. We begin with the description of the opening, the main principle of peak elimination, and its pseudocode. Later, we enrich this algorithm by feature of computing the size spectrum in a single image scan. We also present the arbitrary angle orientation on 2-D image support. The proposed algorithm is based on usage of a queue, which is a FIFO (First In, First out) memory structure. In addition to the basic FIFO features push(), pop(), queue provides front() and back() operations to access the oldest and most recent values, respectively. Each element stored in the queue is composed of two attributes {F, rp}: the pixel graytone value F = f (rp) and its reading position rp in input data stream. Both attributes can be accessed separately, i.e., F x = Q.back() [START_REF] Matheron | Random sets and integral geometry[END_REF] and rp x = Q.back() [START_REF] Maragos | Pattern spectrum and multiscale shape representation[END_REF]. The behavior of this algorithm stems from the principle of a peak elimination. A point of an input signal f (x) is a peak if both its very precedent f (x-1) and subsequent f (x + 1) points are smaller [START_REF] Urbach | Efficient 2-d grayscale morphological transformations with arbitrary flat structuring elements[END_REF]. Our algorithm recognizes 4 possible configurations of these three points ((a) and (b) are peaks, (c) and (d) are not), see Fig. 1, each of which is treated in a different manner.

t f(t) x+1 x x-1 t f(t) x+1 x x-1 (a) (b) t f(t) x+1 x x-1 t f(t) x+1 x x-1 (c) (d)
f (x) > f (x -1) and f (x) > f (x + 1) (7)
The proposed algorithm, see Alg. 1 for the pseudocode listing, computes the open-

ing y = γ B f , y, f ⊂ Z 2 → R, B ⊂ Z is a line segment of length L. It reads input pixel F = f (rp)
and outputs pixel Y = y(rp -L) of the opened image at the time. The points f (x -1), f (x), and f (x + 1), which are needed for the peak elimination, are reachable in the queue as follows: [START_REF] Matheron | Random sets and integral geometry[END_REF] (second to the most recent)

• f (x -1) → Q.back(2)
• f (x) → Q.back(1)[1] (most recent) • f (x + 1) → F = f (rp) (current)
The peak elimination step proceeds in one main while loop (code line 1) that ensures the condition f (x + 1) < f (x) of [START_REF] Urbach | Efficient 2-d grayscale morphological transformations with arbitrary flat structuring elements[END_REF]. If two consecutive pixels are equal f (x + 1) = f (x), the first one is erased from the queue (line 2-4), and replaced by the second later (on line 12). As a consequence, a flat plateau (zone of constant value) is represented in the queue by the last pixel and its position.

Then the second condition f (x -1) > f (x) of ( 7) is tested (Q.back(2)[1]< Q.back(1) [START_REF] Matheron | Random sets and integral geometry[END_REF] on line 6). If the result is false, the f (x) is not a peak and the elimination loop is quit (configuration (d), line 14). Otherwise, the f (x) is a peak and will then accumulated in a single variable. Then the size spectrum S is an array, such as

S[1..L -1] → R +
The size is determined by the distance between f (x) and f (x -1) that need not necessarily be 1 but varies from 1 to L-1 (f (x ± 1) designates precedent/subsequent value of f (x)). The height of a peak obviously depends on the mutual position of f (x + 1) and f (x -1) according to conf. (a), (b), see Fig. 2. The increment of an appropriate size in the size spectrum is carried out on lines 8 and 11.

Arbitrary SE Orientation

The described algorithm can be used for a 2-D input image support as well. It only needs an image to be "mapped" into independent, 1 pixel thin discrete lines (called corridors) oriented in the same angle as the SE. See Fig. 3 for examples of such mapping. Each corridor is then computed according to 1D OPEN in its own queue memory. The major advantage of mapping an image into corridors is that every pixel of input image is read just once, and the result pixel is written within a fixed delay. Furthermore, the input and output data are ordered in the horizontal scan pattern. It allows a stream, in-place, and, using a dedicated hardware, even parallel processing. The corridors are mapped on image by either of two cases in dependence on α. For α ≥ 45 • and α ≤ 135 • , pixel [j,i] (i.e., at line j, column i) belongs to corridor Q ij = (ij sin α)modulo(N + L sin α), where M + L sin α defines the number of queues. The current reading position rp within a corridor is equal to the line index j. On the other hand, for α < 45 • and α > 135 • , pixel [j,i] (i.e., at line j, column i) belongs to corridor Q ij = (ji sin α)modulo(N tan α + L cos α), where (N + L) tan α + 1 defines the number of queues, and the column index i denotes the reading position rp.

1 1 M N 1 1 M N 1 1 M N 1 1 M N (a) (b) (c) (d)
Since each corridor has an independent queue, the image can be read in horizontal scan order regardless its orientation. The image is processed in an ordinary horizontal double-loop, such as for j in 1 to M + L sin α iterate for i in 1 to

N + L cos α that calls y(rp -L) ← 1D OPEN(f (rp), rp, L, Q ij ) of Alg. 1.
Let us observe different image configurations with respect to the SE orientation in Fig. 4. In the horizontal case, the scan order and the SE are parallel so the input is fetched in the right order as the algorithm needs. Therefore, a single queue is used. It reads a value at [k,j] and outputs in the same line at [k,i].

In the vertical or inclined case, the unlike scan and SE orientations require multiple queues. For instance, the inclined orientation in Fig. 4 reads a value at position [l,j] and returns a result at [k,i]. All pixels between these two points are stored in the queues implying the computation latency. Nevertheless, it is the minimal achievable latency considering unlike orientations of the SE and the scan order. The proposed algorithm has very limited memory consumption. The only memory elements are queues whose depth and count are inferred by the SE orientation and image width N . For example, let us consider an opening of 8-bit SVGA image, i.e., 800×600 px by a SE of 41 px. The queues occupy only 656 bits for horizontal, and 525 kbits for vertical SE, respectively, compared to the mere size of the image 3.66 Mbits which is never stored.

Experimental Results

We present the execution time measures to illustrate the computational complexity in this section. We compare the proposed algorithm with a few other efficient algorithms, namely Soille et al. [START_REF] Soille | Recursive implementation of erosions and dilations along discrete lines at arbitrary angles[END_REF], Van Droogenbroeck [START_REF] Van Droogenbroeck | Morphological erosions and openings: Fast algorithms based on anchors[END_REF], and Urbach and Wilkinson [START_REF] Urbach | Efficient 2-d grayscale morphological transformations with arbitrary flat structuring elements[END_REF]. The benchmarks were performed on Intel Xeon E5620 @2.4GHz CPU, running Linux. The time reported in tables below refers to the user CPU time consumed by the respective algorithms. We use the mountain natural photo as a testing image, originally introduced at [START_REF] Van Droogenbroeck | Morphological erosions and openings: Fast algorithms based on anchors[END_REF]. At first, we evaluated the execution time benchmark with respect to the increasing image size, see Fig. 5. The results illustrates that the complexity of our algorithm is linear with the image size.

The second benchmark in Fig. 6 retains the same image size and changes the length of the horizontal SE to show that the execution time is independent of SE size (however the value of execution time is dependent on the image content). Compared with other algorithms, only Van Droogenbroeck's opening algorithm outperforms our solution. On the other hand, our algorithm computes the size spectrum in addition to the opening and is independent of the number of gray levels. For example, let us consider the size spectrum up to λ MAX = 100. Even if we omit the arithmetic operations of (4), the pure time for computation λ MAX γ λ will take 100×2.7 ms = 270 ms using the fastest Van Droogenbroeck's algorithm. Our algorithm computes the size spectrum in a single run, i.e., in 9.9 ms.

Applications

This section deals with examples of suitable applications of arbitrary oriented opening and size spectrum. The first one is the local orientation information (orientation field) ζ 35 of a fingerprint image in Fig. 8. The same operator can be applied to detect angle of the road lines in Fig. 8, too. Later, Fig. 9 displays the oriented size spectra OS(α, λ) for two different materials.

Conclusions

This paper presents a new algorithm for 1-D morphological opening and size spectrum extraction with large scale of suitable applications, such as anisotropic texture analysis, size distribution, local orientation, image restoration, etc.. The algorithm processes an image in a constant time regardless the SE size and almost regardless its rotation. The operator always uses strictly sequential access to data in horizontal scan pattern. Furthermore, it infers the minimal latency and very low memory requirements. The proposed algorithm achieves a high performance in the sequential C implementation and outperforms other arbitrary oriented line opening algorithms. Moreover, it computes the size spectrum in a single image scan introducing a speed-up in orders of magnitude compared to residue approach demanding multiple opening runs. The future extension of this work is to implement the proposed algorithm in the FPGA circuit taking an advantage of both spatial and temporal parallelism.

Figure 1 :

 1 Figure 1: Four different pixel configuration for peak identification. Configuration (a) and (b) characterize a peak, whereas configuration (c) and (d) do not.

Figure 2 :

 2 Figure 2: Size spectrum increment for (a) configuration a, and (b) configuration b. (c) The peak is sliced by each gray level it contains.

Figure 3 :

 3 Figure 3: Image corridors (discrete lines) mapping for different SE orientations (a) a horizontal SE, (b) a vertical SE, and (c) and (d) inclined SEs.

Figure 4 :

 4 Figure 4: Image configuration for different SE orientations (a) a horizontal SE, (b) a vertical SE, and (c) an inclined SE. × denotes the pixels to be read in next iterations, • denotes the previous output pixels.

→

  Execution time [s]

Figure 5 :

 5 Figure 5: Execution time of opening/spectrum versus the image size. Structuring element is vertical 101 px.

→

  Size of horizontal SE [px]→ Execution time[s]

Figure 6 :→

 6 Figure 6: Execution time of opening versus the size of the horizontal structuring element. Natural photo 800 × 600 px is used.

Figure 7 :

 7 Figure 7: Execution time of opening versus the rotation angle α. Natural photo 800 × 600 px is used.

Figure 8 :

 8 Figure 8: Extraction of local orientation on (a) a fingerprint and (b) a road image.

Figure 9 :

 9 Figure 9: Oriented size spectrum OS(α, λ). △ α = 1 • , λ MAX = {50, 100}.

). This is decided upon condition f (x + 1) < f (x -1) (F < Q.back(2) [START_REF] Matheron | Random sets and integral geometry[END_REF] on line 3). Obviously, the while loop iterates until a non-peak configuration ((c) or (d)) is encountered. When all peaks are erased, the current pixel value is unconditionally pushed into the queue along with current reading position (line 15). The oldest stored pixel is checked whether it has been stored in the queue for too long. This check is carried out by comparing the stored reading position plus SE size with current rp (line 16). Outdated values are immediately deleted. The oldest stored value Q.front(1) [START_REF] Matheron | Random sets and integral geometry[END_REF] is the result of Alg. 1 as soon as rp exceeds the SE size (line 18).

The algorithm presented so far computes the opening transformation using the principle of the peak elimination that discards the peak values. When we extend the algorithm by a feature for storing the eliminated peaks, we obtain the size spectrum S with a minimal additional effort. As the algorithm eliminates a peak recursively, it literally slices the peak by each gray level of a peak. The slices of the same size are