Ahmed Zermane

Niels Dehio

Fellow IEEE Abderrahmane Kheddar

Planning Impact-Driven Logistic Tasks

Keywords: Planning Impact-Driven Logistic Tasks Logistics, Impact-aware planning and control, Tossing, Grabbing, Boxing, BiRRT

HAL is

I. INTRODUCTION

I MPACT-DRIVEN tasks are challenging to plan and to control due to non-smooth and fast contact dynamics. We are particularly interested in tasks found in the automated distribution and sorting logistics chains. Namely, fast and swift picking and tossing, grabbing, and boxing of objects of all sorts, see Fig. 1. The challenge is that the manipulated objects can undergo (open-loop) phases that are out of the robot control. For example, tossing an object requires accurate computation of its ballistic trajectory to meet landing specifications. Another example is manipulating (boxing) an object with a robot using impulses inducing an initial velocity (or, equivalently, an energy input) to bring it from an initial pose to a desired one: after the impulse, the object's motion is not 'controlled' by the robot. Our main contribution consists in devising a strategy in terms of planning and control of such tasks (defined by project partner Vanderlande) in two separate modules (Fig. 1):

• M1: Plan for the robots a trajectory (with or without a load) to reach a target location (position and orientation) and velocity (linear and angular) from a given initial state (comprising both position and velocity). We allow the starting velocity to be non-zero for the robot not need to stop at each cycle. Our planning module accounts for hardware and taskspace limits; see Sec. III.

• M2: For each task (tossing in Sec. IV, grabbing in Sec. V, Fig. 1: Proposed motion generator for scenarios with impact. and boxing in Sec VI), we devise separate modular plugins representing specific inverse kinematics solvers for mapping desired impact objectives from task-space onto joint-space. The output of M2 is the input of M1, i.e., the planner ensures bypassing the via-point with desired position and velocity.

With respect to related work (discussed in Sec. II), our novelty lies in defining task coordinates as a via-point, which are integrated into a time-constrained path. The latter is optimized for a generated samples' map. Throughout the entire letter, we assume rigid objects with known geometry and wellestimated inertial parameters [START_REF] Kurdas | Online payload identification for tactile robots using the momentum observer[END_REF].

We assessed our approach in several real-robot experiments, finally demonstrating a complex industry-driven task that could be customized to any sequence order, see Sec. VII.

II. A BRIEF BACKGROUND

• Planning and motion generation are well-explored topics in robotics. A multi-mode trajectory optimization proposed in [START_REF] Stouraitis | Multi-mode trajectory optimization for impact-aware manipulation[END_REF] combines hybrid control and dynamics to embrace impact for stopping moving objects in impact-aware manipulation tasks. In [START_REF] Toussaint | Describing physics for physical reasoning: Force-based sequential manipulation planning[END_REF], a multi-physics formulation enables one to interact with objects at non-zero relative velocity to achieve desired motions. Alternatively, by relying on trajectory optimization under contact constraints, the robot can push objects [START_REF] Sleiman | Contact-implicit trajectory optimization for dynamic object manipulation[END_REF]. Whereas [START_REF] Mészáros | Learning to pick at non-zerovelocity from interactive demonstrations[END_REF] relies on learning from interactive demonstrations to execute picking at non-zero relative velocities. In [START_REF] Wang | Impact-aware taskspace quadratic-programming control[END_REF], an impact-aware QP controller is formulated to be resilient to 3D frictional impacts in the task-space by defining new impact-related constraints and reformulating the usual ones to turn them robust against abrupt changes of the states at impacts. Other related works extend to impact for different robots and tasks within [START_REF] Khurana | Learning to hit: A statistical dynamical system based approach[END_REF], [START_REF] Bombile | Dual-arm control for coordinated fast grabbing and tossing of an object: Proposing a new approach[END_REF] using dynamical systems and machine learning-based approaches.

• Tossing can be seen as an extension of the robot's workspace or for speeding pick-and-place task that do not require high placement precisions. Existing robotic tossing methods are found in model-based and data-driven approaches. The most outstanding work is reported in [START_REF] Zeng | Tossing-Bot: Learning to throw arbitrary objects with residual physics[END_REF]. It is a machine learning-based system combined with residual physics designed to enable rapid picks and throws of arbitrary objects into selected boxes. In [START_REF] Kasaei | Throwing objects into a moving basket while avoiding obstacles[END_REF], tossing is formulated using reinforcement learning (RL), enabling generalization on various objects and fast reaction to dynamic environments. Combining decision transformers (DT) and RL enables tossing objects from handful experiments and far-from-reality trials [START_REF] Monastirsky | Learning to throw with a handful of samples using decision transformers[END_REF]. Earlier work in [START_REF] Ghadirzadeh | Deep predictive policy training using reinforcement learning[END_REF] handled basic shapes (balls) grasping and tossing using a deep neural network architecture to perform grasping. Throwing a ball by a kinetic chain approach based on swinging motion and constraining its motion to the release timing and release direction is proposed in [START_REF] Senoo | High-speed throwing motion based on kinetic chain approach[END_REF]. A dedicated gripper for grasping and throwing is proposed in [START_REF] Govindan | A novel hybrid gripper capable of grasping and throwing manipulation[END_REF], containing a latching mechanism that ensures both high accelerations for throwing by instantly discharging the stored elastic potential and low for placing.

• Dual-am grabbing enables higher payload handling and larger object grasp, e.g., [START_REF] Smith | Dual arm manipulation-a survey[END_REF]. In a previous work [START_REF] Dehio | Dual-arm box grabbing with impact-aware mpc utilizing soft deformable end-effector pads[END_REF], this scenario was achieved using an impact-aware model preview control with learned dynamics of soft deformable pads mounted on the end-effectors. In [START_REF] Bombile | Dual-arm control for coordinated fast grabbing and tossing of an object: Proposing a new approach[END_REF], Modulated Dynamical Systems (MDS) generate motion for a dual-arm system to grab and release an object reactively swiftly. In [START_REF] Van Steen | Dual arm impact-aware grasping through time-invariant reference spreading control[END_REF], a time-invariant reference spreading controller is formulated to manipulate the object swiftly. Impact grabbing is studied more fundamentally and theoretically in [START_REF] Aydinoglu | Stabilization of complementarity systems via contact-aware controllers[END_REF], [START_REF] Brogliato | Modeling, analysis and control of robot-object nonsmooth underactuated lagrangian systems: A tutorial overview and perspectives[END_REF], where Linear Complementarity Systems (LCS) formulation is employed to synthesize controllers and study systems' stability in rich contact manipulation problems, and dual-grabbing as an example. In [START_REF] Shirai | Chanceconstrained optimization in contact-rich systems[END_REF], grabbing is treated as a Stochastic Discrete-time LCS (SDLCS) solved by a chance-constrained based optimization.

• Boxing, i.e., pushing and tilting objects can be treated as impact-contacts in motion control and planning [START_REF] Mason | Toward robotic manipulation[END_REF]. For example, [START_REF] Cheng | Manipulation with suction cups using external contacts[END_REF] used a suction cup end-effector to tilt an object by developing a framework that withstands deformation modeling errors. In [START_REF] Shirai | Robust pivoting manipulation using contact implicit bilevel optimization[END_REF], a robust pivoting formalism accounts for object inertia parameters. It is developed based on contact implicit bilevel optimization and enhanced by a tactile closedloop control to deal with initializing issues and enable failure recovery. Impact-based interactions based on contact-implicit trajectory optimization (TO) are proposed in [START_REF] Sleiman | Contact-implicit trajectory optimization for dynamic object manipulation[END_REF], assuming a hard contact model and imposing unilateral constraints as complementarity conditions. On the other hand [START_REF] Onol | Contact-implicit trajectory optimization based on a variable smooth contact model and successive convexification[END_REF] proposed a TO based on a variable smooth contact model and successive convexification to push objects to a desired pose, whereas for the same problem [START_REF] Khurana | Learning to hit: A statistical dynamical system based approach[END_REF] used a learned statistical dynamical system to generate a non-zero relative contact velocity.

III. IMPACT-DRIVEN TRAJECTORY PLANNING

Module M1 builds on a sampling-based algorithm Bi-RTT [START_REF] Kuffner | Rrt-connect: An efficient approach to singlequery path planning[END_REF]. It uses a steering function to plan minimum-time trajectories in the robot's configuration space subject to a viapoint constraints in the task-space. The latter embeds explicitly impact postures and enforces back-to-safety configurations. Unlike any of the related works (e.g., [START_REF] Berscheid | Jerk-limited real-time trajectory generation with arbitrary target states[END_REF]), this creates profile motions with way-points on position, velocity and acceleration levels. This Module relies on acceleration-limited trajectory generation where Double Integrator Minimum Time (DIMT) is the neighborhood criterion [START_REF] Kunz | Probabilistically complete kinodynamic planning for robot manipulators with acceleration limits[END_REF]. It allows having non-zero start and end velocities of the planned trajectories. If any jerk constraints exist, it can be used within the steering function. Fig. 2: Two grown trees (blue and red edges). Sampled states are chosen within the bounds for each DoF, defined by velocity, position, and acceleration limits in the jerk-limited case. Via-point constraints for a task-space pose and velocity are converted into joint-space states using inverse kinematics.

To reach a target position/velocity (pose) state from a given initial one, Bi-RTT grows iteratively two trees: (i) forward in time starting from the initial state, and (ii) backward in time starting from the final state, see Fig. 2, blue-dotted area on the planning flowchart in Fig. 3. The sample-map segments with limited duration prevent numerical integration errors. Hence, there is no need to define an upper limit time [START_REF] Berscheid | Jerk-limited real-time trajectory generation with arbitrary target states[END_REF] for the whole trajectory. For the acceleration-limited trajectory generation, each state Q is represented by the stacked vectors of joint states Q = [q, q], which extend to accelerations Q = [q, q, q] in the jerk-limited case.

The steering method, described in Algorithm 1 (lines 12 and 18) and detailed in Algorithm 2, is either the DIMT or the jerk-based one, see details later. This steering method relies on determining the closest neighbor using a distance function. Different metrics exist to define distance functions. We use the go-to time between two states [START_REF] Kunz | Probabilistically complete kinodynamic planning for robot manipulators with acceleration limits[END_REF], [START_REF] Hauser | Fast smoothing of manipulator trajectories using optimal bounded-acceleration shortcuts[END_REF]. Despite not being a metric because of its asymmetry, it is proven reliable, parameter-free, and accounts for both position and velocity. As shown in Algorithm 1, a connection between two trees can be made if the randomly generated state Q r is successfully added to both trees (lines [START_REF] Dehio | Dual-arm box grabbing with impact-aware mpc utilizing soft deformable end-effector pads[END_REF][START_REF] Van Steen | Dual arm impact-aware grasping through time-invariant reference spreading control[END_REF][START_REF] Aydinoglu | Stabilization of complementarity systems via contact-aware controllers[END_REF][START_REF] Brogliato | Modeling, analysis and control of robot-object nonsmooth underactuated lagrangian systems: A tutorial overview and perspectives[END_REF][START_REF] Shirai | Chanceconstrained optimization in contact-rich systems[END_REF][START_REF] Mason | Toward robotic manipulation[END_REF][START_REF] Cheng | Manipulation with suction cups using external contacts[END_REF][START_REF] Shirai | Robust pivoting manipulation using contact implicit bilevel optimization[END_REF].

Each trajectory task can be generated by skipping path selection at the first phase of sample-map building for N iterations. For the same task, each via-point has multiple jointspace configurations competing to generate the best path, reddotted area in the flowchart see Fig. 3. We start evaluating the desired number of IK solutions (loop-line 9) by considering it as the new random states (line 10). Finally, Algorithm 1 extracts the path that fulfills all the initial and final conditions under limits and collision constraints (line 27). It has a time complexity of O(2n), where n = max(N (V 1 , V 2 )). To introduce both steering functions, we define:

Algorithm 1 Extended BiRRT Algorithm 1: procedure BIRRT(qmax, q min , qmax, qmax, P init , P final , N) 2: cost i ← 0, cost f ← 0 3: FinalPath ← [ ]; FinalCost ← ∞; 4: V 1 ← [P init , Q init , cost i ]; 5: V 2 ← [P final , Q final , cost f ]; 6: E1 ← [ ]; E2 ← [ ]; 7:
d ← true; % growing direction indicator 8:

counter ← 0; 9:

while counter ≤ N do 10:

Qr ← SampleInJointSpace; 11:

Pr ← ForwardKinematics(Qr), CollisionCheck(Qr); 12:

T 1 , S 1 , costr i , JointLimits ← Steering(Qr, V 1, d); 13:
if ¬JointLimits then 14:

Q i , P i ← SegmentsGeneration(T 1 , S 1 ); 15: Etemp ← BuildEdges(Q i , P i , S 1 ); 16: V 1 ← V 1 ∪ [Pr, Qr, costr] ∪ [P i , Q i , cost i ]; 17: E1 ← E1 ∪ Etemp; 18: T 2 , S 2 , costr f , JointLimits ← steering(Qr, V 2, d); 19:
if ¬JointLimits then 20: • ( qmax , qmax , ... q max ) velocity, acceleration, jerk bounds.

Q i , P i ← SegementsGeneration(T 2 , S 2 ); 21: Etemp ← BuildEdges(Q i , P i , S 2 ) ; 22: V 2 ← V 2 ∪ [Pr, Qr, costr] ∪ [P i , Q i , cost i ]; 23: E2 ← E2 ∪ Etemp; 24: Cost Tmp = costr i + costr f ; 25: if Cost Tmp < FinalCost then 26: FinalCost ← Cost Tmp ; 27: FinalPath ← Extractpath(V 1 , E 1 , V 2 , E 2 , d); 28: end if 29: end if 30: end if 31: Swap([V 1 , E 1 ], [V 2 , E 2
In DIMT steering function, there are four possible extremal motions for velocity and acceleration bounds, see Fig. 4: the parabolas P + and P -accelerating at qmax and at -q max , respectively; and the lines L + and L -traveling with zero acceleration at max velocity qmax and at min velocity -qmax , respectively. The sign of the first acceleration segment is determined programmatically [START_REF] Kunz | Probabilistically complete kinodynamic planning for robot manipulators with acceleration limits[END_REF], using σ = sgn(q i+1 -q i -δq acc ),

where δq acc is the distance crossed while accelerating at qmax , which is defined as: Fig. 4: Position-velocity-Acceleration profiles, Two segments profile, and Three segments profile.

δq acc = 1 2 ( qi + qi+1 ) | qi -qi+1 | qmax . (2)
thus we have qi1 = -q i2 = σ qmax and qlimit = σ qmax . The motion profile time parameters are defined as follows: 1) determine t 1 by solving the quadratic equation,

qi1 t 2 1 + 2 qi t 1 + q2 i+1 -q2 i 2q i2 -(q i+1 -q i ) = 0; (3)
2) check for velocity bounds violation, if ( qt1 < qmin ) or ( qt1 > qmax ), then it is a trapezoidal constantvelocity profile segment with duration t v :

t 1 = qlimit -qi qi1 , (4)
t v = q2 i+1 + q2 i -2 q2 limit 2 qlimit qi1 + q i+1 -q i qlimit , (5)
t 2 = qi -qlimit qi2 . (6)
else:

t 2 = qi+1 -qi qi2 + t 1 . (7)
After obtaining the extremal profile for each DoF with duration T i (i = 1, n and n is DoFs number), and checking for its blocked time interval, the minimum-acceleration trajectory [START_REF] Kunz | Probabilistically complete kinodynamic planning for robot manipulators with acceleration limits[END_REF], [START_REF] Lau | Smooth RRT-connect: An extension of RRTconnect for practical use in robots[END_REF] is then applied to synchronize T for all the DoFs by calculating new profiles parameters as follows:

T 2 qi1 + (2T ( qi+1 + qi ) -4(q i+1 -q i ))q i1 -( qi+1 -qi ) 2 = 0, (8)
For the Jerk-limited steering method. We use the algorithm proposed in [START_REF] Berscheid | Jerk-limited real-time trajectory generation with arbitrary target states[END_REF] (available online). The resultant motion profile is an S-curve time scaling that consists of a maximum of seven segments.

IV. TOSSING Our approach differs from all previously cited work [START_REF] Zeng | Tossing-Bot: Learning to throw arbitrary objects with residual physics[END_REF], [START_REF] Kasaei | Throwing objects into a moving basket while avoiding obstacles[END_REF], [START_REF] Monastirsky | Learning to throw with a handful of samples using decision transformers[END_REF], [START_REF] Senoo | High-speed throwing motion based on kinetic chain approach[END_REF], [START_REF] Govindan | A novel hybrid gripper capable of grasping and throwing manipulation[END_REF] in explicitly defining constraints reflecting any robot's ability to toss within a computed reach workspace, and more importantly, in integrating explicit constraints on objects' impact specifications. That is to say, we account for objects landing impact location, relative speed, or impulse. Then, it is plugged into the developed planner previously with the physical modeling of the tossing ballistic.

Prior to planning and control, we define the toss release reachable workspace (in both position and velocity). The Algorithm 3 describes the entire process, and Fig. 5 illustrates the outcome in a schematic way.

Algorithm 3 Tossing-Workspace (Robot, Z max , Z min , step) 1: RWspace ← ReachableWorkspace(Robot, V, d); 2: z = Z min ; 3: while z ≤ Zmax do 4:
TWspace ← TossingWorkspace(Robot, RWspace, z); 5: z = z + step; 6: end while 7: return TWspace; First, we define the robot's reachable workspace with task velocity bounds (line 1 Algorithm 3). Using the latter, we estimate the tossing workspace by solving simple ballistic motion discarding post-impact effects (e.g., bouncing, sliding...) for each height layer (z). For each z we sample a launching position with different velocities (using an iterative algorithm based on Monte Carlo). We solve for each pose (position and velocity) (line 4 Algorithm 3) the landing time t l such that P t l = z l from eq. ( 9) (see later). Fig. 5: example of a reachable workspace (sphere-like volume with black edges) and tossing workspace, for z max = 0.5 m and z min = -0.7 m (left). Tossing set for a given target (right).

We formulate the tossing ballistics as an optimization problem. The release state (a via-point in our planning) is determined by optimizing the ballistic motion, considering the robot's capabilities to generate a feasible one. Other ballistic related phenomena could be modeled, e.g., aerodynamics, bouncing, gripper dynamics... but would lead to a complex non-linear problem. Our primary goal is to solve tossing from industry-driven perspectives. Bouncing is challenging to predict. Therefore, one of the objectives we set is having low impact relative-velocity of the thrown objects when the precision of the landing spots prevails.

We define

P d ∈ R 3 , Q d ∈ R 4
as the desired position and quaternion orientation of the object and P t l , Q t l the ones at time-step t l of the ballistic motion defined by:

P t l =   x y z   =   v x t l + x o v y t l + y o -gt 2 l + v z t l + z o   , (9)
and

Q t l = t l 0 1 2 Ω(ω)Q dt with ω = I -1 τ -I -1 (ω × Iω) (10)
where Ω(ω) ∈ R 4×4 is a skew-matrix; g is the gravitational acceleration; [x o , y o , z o ] ∈ R 3 and Q o ∈ R 4 are position and orientation coordinates of the object's CoM at the moment of release; v ee = [v x , v y , v z ], ω ee = [ω xo , ω yo , ω zo ] are the linear and angular velocities of the end-effector from which that of the picked object is determined as initial conditions for eq. ( 9) and eq. ( 10), respectively; I ∈ R 3×3 is the object inertia matrix; and τ ∈ R 3 is the external torques for a free ballistic motion (assumed nil, τ = 0 3×1 ). We denote the translation error e P = P t l -P d ; used as design, evaluation criterion in [START_REF] Zeng | Tossing-Bot: Learning to throw arbitrary objects with residual physics[END_REF], [START_REF] Monastirsky | Learning to throw with a handful of samples using decision transformers[END_REF] respectively; and express the orientation error e O as the imaginary part of the Hamilton product between desired and landing rotations:

e O = [Q d Q -1 t l ] = ν t l • d -ν d • f l + d × t f ( 11
)
The problem of the IK for tossing accounting for impact constraints (location, orientation, and velocity) and robot limits, is formulated as a weighted optimization as follows: min q, q w 1 e P 2 + w 2 e O 2 + w 3 v desired -v impact 2 s.t: qmin ≤ q ≤ qmax and q min ≤ q ≤ q max [v ee , ω ee ] = J(q) q v ee ≤ v max and ω ee ≤ ω max [START_REF] Ghadirzadeh | Deep predictive policy training using reinforcement learning[END_REF] where decision variables (q, q) are robot's joint position and velocity; J is the Jacobian matrix; and w i the associated weight to each term. The object's impact velocity v impact is considered in the ballistic optimization. The generated tossing pose, chosen from a constrained subset in Fig. 5, would induce a landing velocity as close as possible to v desired at the landing time t l , i.e., by the end of the ballistic motion trajectory. The solution of the optimization problem defined in eq. 12 is then considered a via-point state for the trajectory planner. The target {P d , Q d , v impact } belongs to the tossing workspace; else the solution found will be the closest to it.

V. IMPACT GRABBING Swift grabbing of objects with impact using a dual-arm robotic system requires (i) knowing where to contact the object for each arm (e.g., could be driven by having the object balanced during lifting), (ii) specifying relative velocity at both impacts, and, more importantly, (iii) synchronizing both arms to impacts at the same time. To fulfill these requirements, we query the planner to find feasible collision-free configuration paths for both robots (which can be extended to multi-robots). Indeed, our proposed planner can straightforwardly handle dual manipulation and many configurations can be put into cascade to form a fully defined series of tasks, i.e., dual-grab, dual-toss, or dual-place, dual-re-grab, etc.

Synchronizing both robots can be done at the planning level; time constraints are added to the planned path for each robot to match the arrival time for both of them. It can also be made at the control level since we use the same tracking strategy in Sec. VII-A, introducing the posture task of the fastest robot with a delay corresponding to the time difference between the planned paths.

Once the impact with the object occurs, the robot shall instantly switch the control objectives in terms of forces and torques toward a stable grasp and swift motion control. Our multi-robot QP task-space control formulation enables parallel admittance tasks to monitor and control the contact state between the dual-arm system and the object, see Fig. 6. In order to account for impact uncertainties due to one robot contacting the object before the other one, we mitigate the impact velocity bounds (using momentum conservation) in the impact-aware QP controller w.r.t to the possible motion of the object in the direction of the opposite contact, which results in lowering impacting velocity bounds for both robots.

VI. IMPACT BOXING

Boxing involves applying an impulse on an object to induce an initial velocity that brings it to a desired state. Concretely, we examine how to rotate a resting box using a single impulse provided by a robotic arm, which is mainly different from what is described in related work, e.g., [START_REF] Sleiman | Contact-implicit trajectory optimization for dynamic object manipulation[END_REF], [START_REF] Khurana | Learning to hit: A statistical dynamical system based approach[END_REF], [START_REF] Mason | Toward robotic manipulation[END_REF], [START_REF] Cheng | Manipulation with suction cups using external contacts[END_REF], [START_REF] Shirai | Robust pivoting manipulation using contact implicit bilevel optimization[END_REF], [START_REF] Onol | Contact-implicit trajectory optimization based on a variable smooth contact model and successive convexification[END_REF]. We use an explicit model-based plugin that integrates our planning strategy by accounting for all forces acting on the box; see Fig. At impact time T , utilizing Newton's Laws circumvents the complexity linked to state discontinuities [START_REF] Halm | A quasi-static model and simulation approach for pushing, grasping, and jamming[END_REF]. Thus, there is a need to use second-order Measure Differential Equation (MDE) (assuming no sliding):

I yy /o 1 ω = τ ι δ T j + τ P j (13)
and

τ ι = # ‰ OM × F ι = [0, -z M f ιx + x M f ιz , 0] T τ P = # ‰ OG × P = [0, mgx G , 0] T = [0, mgd 1 cos(α 1 ), 0] # ‰ OM = l cot(ψ), 0,l T , # ‰ OG = d 1 cos(α 1 ), 0, d 1 sin(α 1 )
are the coordinates of the impact point M and the object's CoM; F ι is the impact force generated by the robot endeffector; P is the gravity force acting on the box; and δ T is the Dirac measure at impact. [I, ω] is described in Sec. IV. We enforce non-sliding impact with

| f t | ≤ µ 0 f n , | τ ι | ≥ D y f n , andf ιz > f ιx > 0 ( 14
)
where µ 0 and D y are friction coefficients. Post-impact dynamics can be expressed as

I yy /o 1 ω = τ P (θ) j = mgd 1 cos(θ) j ( 15
)
with ω 0 = ω + 0 = 0 resulting from the impact, and θ + = 0. the post-impact velocity at M is written as: [START_REF] Dehio | Dual-arm box grabbing with impact-aware mpc utilizing soft deformable end-effector pads[END_REF] we get ω + by solving the energy based equation eq. 17:

ν + M = ω + × # ‰ OM = ω + [l cot(ψ), 0, l] T (
1 2 I yy /o 1 ω + 2 = mgd 1 (1 -sin(α 1 )) (17
)

ω + = 2mgd 1 I yy /o 1 (1 -sin(α 1 )) ( 18
)
we can then estimate the impact point and velocity from eq. 16 to plug it into the impact planner. However there are situations where, even in perfect physics, it is not possible to reach a settled pose. The right amount of velocity to bring the object to reach its intended pose induces a rotational trajectory from the initial resting pose to the singular point at near-to-nil angular velocity (S 1 in Fig. 7). Having the right minimal energy making the rotating object pass the singular pose S 1 , the remaining second phase trajectory is a free-falling one; which induces another impact at S 2 resulting in an energy gain. Depending on the impact physics, the geometry of the object and the landing ground, the desired pose can be reached after a few bounces or the object would fall on the other side.

In order to achieve boxing, some conditions must be satisfied, first

τ ι > τ P ⇔ l -f ιx + cot(ψ)f ιz > mgd 1 cos(α 1 ) (19)
From eqs. ( 14) and ( 19), we constraint x M

f ιx < cot(ψ)f ιz ⇔ ψ > π 4 ⇔ x m < l (20)
Then from eq. 19, the influence of the box thickness is also clear, i.e., the thinner the box to be flipped, the more impact force is needed. At some point it will never be possible to flip it with any amount of normal force

l → 0 ⇒ x m → 0 ⇒ ψ → π 2
The CoM position to the rotation axis can also make this task difficult to achieve. By analyzing energy at singularity positions S 1 and S 2 (see Fig. 7), we seet two main conditions:

(i) θ 2 + α 2 ≤ π 2 , and (ii) E 2 = mgd 2 cos(θ 2 + α 2 ) ≤ E 1 = mgd 1 cos(α 1 )
, where is the coefficient of restitution.

Assuming that, at S 1 , ω S1 ≥ 0; at S 2 , ω S2 = 0; and = 1,

mgd 1 = mgd 2 cos(θ 2 + α 2 ) (21)
since cos(

θ 2 + α 2 ) ∈ [0, 1] we have d 1 ≤ d 2 (22)
which means π 2 -α 2 ≤ α 1 .

VII. EXPERIMENTS

We conducted extensive experiments using 7 DOF Panda robots to assess our approach. We also used the open-source QP controller mc_rtc that runs at 1 kHz using the necessary dependencies [START_REF] Dehio | Dual-arm box grabbing with impact-aware mpc utilizing soft deformable end-effector pads[END_REF] mc_franka and mc_panda. We report some results and data obtained and provide an accompanying video Check video: https://youtu.be/A8N0C3mE3TM as supplementation material.

A. Trajectory Tracking

Common to all tasks, the planner (Sec. VII-A) does not include the robot's dynamics. Hence, we cannot track the planned trajectory obtained from Algorithm 1 in open-loop; therefore, we impose a posture task objective in the QP controller [START_REF] Bouyarmane | On weight-prioritized multitask control of humanoid robots[END_REF]. Let the triplet {q ref , qref , qref } t be position, velocity, and acceleration references at t time-step within the generated trajectory, and {q cur , qcur } t the current position and velocity sensor readings. The desired task-space acceleration is formulated as a classical linear spring-damper feedback law qd = qref -K s (q cur -q ref ) -K d ( qcur -qref ), with diagonal positive-definite matrices K s , K d that represent the stiffness and damping gains. The QP is formulated as:

min q∈R n q -qd s.t: C( q, q)q + d ≤ 0 (23)
where C( q, q) is the constraints vector related to the robot's limits (position, velocity, acceleration, torque, jerk); the jerk is included in the constraints by deriving generated acceleration; d is the bounds vector. Figure 8 shows a sample trajectory-tracking result, revealing a limited position and velocity error: |e position | ≤ 6.10 -3 rad, |e velocity | ≤ 25.10 -3 rad/s for all joints, and the corresponding task-sapce velocity and position errors to the via-point. Which is enough for the intended tasks and their application in a logistics context under the assumption made for each task.

B. Tossing

First, the robot's reachable and tossing workspaces are constructed. Then assuming that we know precisely the object inertia parameters, we toss to several chosen empty containers as landing targets (put at different heights) within the reachable tossing space. These tossing experiments are achieved either with a suction cup or a panda gripper as picking tools. For the suction cup case, the object is assumed to be rigidly attached to the suction cup for the planning, and the suction cup -object contact rupture delay is determined statistically by analyzing pressure drop down and actual release command at around ≈ 200 ms in preliminary experiments. This is to show that we are able to plan various objects tossing on different locations. We also specified different landing velocities (equivalently different landing impacts) w.r.t to the surface normal for the same landing target of the same object. To highlight differences between ballistic motions, we compare the release configurations and measure the object impacts. For this purpose, we instrumented the target floor with ATI-force sensors; see result data in Fig. 9. The obtained results demonstrate that we are able to toss objects (providing their inertia parameters are known) precisely to any reachable location surrounding the robot. Yet, contrarily to existing tossing approaches, see Sec. II (tossing) and Sec. IV (first paragraph), we are able to put specifications on the landing pose and relative velocity (or impact). We also didn't have to learn from any real or simulation data/trials or put constraints on the tossing configuration (shortcut used in the learning methods) that limits the performance of robots and reduces the solutions set. For now, it is still difficult to precisely predict object pose and location after impact, as such behaviors are yet extremely difficult to predict.

C. Grabbing

We have also replicated intensive grabbing scenarios, namely those already achieved in [START_REF] Dehio | Dual-arm box grabbing with impact-aware mpc utilizing soft deformable end-effector pads[END_REF]. The only addedvalue with respect to our latter work is that now, we are able to automatically generate grabbing trajectories that were made ad-hoc by the user in [START_REF] Dehio | Dual-arm box grabbing with impact-aware mpc utilizing soft deformable end-effector pads[END_REF], see results on grabbing in Sec. VII-E where all skills are combined. Compared to [START_REF] Bombile | Dual-arm control for coordinated fast grabbing and tossing of an object: Proposing a new approach[END_REF], we do not need to have an intermediary MDS learning phase, as we can simply define task-space targets. Moreover, in [START_REF] Bombile | Dual-arm control for coordinated fast grabbing and tossing of an object: Proposing a new approach[END_REF], there is no guarantee on constraint fulfillment; also, it is not possible to specify landing constraints other than object location.

D. Boxing

We tested boxing (tilting fully different boxed-objects) utilizing a 3D-printed rigid end-effector. Boxes are used in our experiments because a large quantity of parcels found in Fig. 11: Desired joint velocity for three joints, i = {2, 3, 4}, dark colors (red, bleu, green) and the real measured corresponding velocities by the encoders with light colors.

Figure 11 illustrates the velocity tracking performance in the boxing experiment (namely in the presence of velocity jumps produced by the impact), and how important it is to take the robot back to a safe configuration to proceed to with the next task objective.

In case of boxes that cannot reach the desired pose in a stand-alone motion (i.e., would bypass S 2 in Fig. 7), it is possible to plan another robot arm to meet the box a little after the singular configuration S 1 in Fig. 7 and still achieve boxing.

There is no similar existing work to which our approach is to be compared. Indeed, the boxing studied in this paper is an extreme dynamic motion case; in most cases, it is acceptable to gently tilt boxes to grab them with slower dynamics by applying a continuously sustained contact. This can be done as formulated in [START_REF] Shirai | Tactile tool manipulation[END_REF] using tactile feedback, where contact is sustained during the whole motion. In terms of performance, rotating a box in [START_REF] Shirai | Tactile tool manipulation[END_REF] took ∼ 160 s. Our implementation of contact-sustained boxing using force sensing took 45 s, whereas impact-boxing is made in around 1 s!

E. Combining different tasks

With our modular impact-aware planning and control approach, it is straightforward to combine previous individual impact-aware tasks. Since the planning can start from any state, no pause is necessary between two successive trajectories. We successfully experimented with a scenario that can be encountered in logistics were a given object at rest would first be boxed to change its pose, then grabbed and then tossed. Figure 12 illustrates such an example with snapshots from the multimedia video of this scenario . For the latter scenario, Fig. 12: Sequencing boxing, grabbing, and tossing in one run. Fig. 13 shows the force measurements for part of the task starting from boxing to dual arm tossing (notice that tossing is made in dual arm, that is adapting the computation of the working spaces) indicating boxing performed by the right panda robot (light red region) and the transition from dual grabbing impact (blue region) to sustained contact through an admittance task objective in the QP controller (orange region). Fig. 13: Measured right panda (see Fig. 12) end-effector applied forces (along x, y, and z-axes) while impacting a box to change its initial rotation. In the grabbing interval, the endeffectors normal forces are shown (also during contact).

F. Discussions and limitations

Although repeatable, our experiments are conducted in wellcalibrated settings, blindly (i.e., without using vision for object tracking). Adding object tracking would enhance significantly the performance and robustness of our approach. Moreover, the shape and inertia parameters of the objects are known, which is not always the case in logistics sorting. Yet, online and fast inertia parameter estimation is possible [START_REF] Kurdas | Online payload identification for tactile robots using the momentum observer[END_REF] and has to be integrated as a prior task in our planning, which is certainly possible for grabbing or tossing. Yet, for boxing, we cannot escape from first impact that would allow identifying such parameters for unknown objects. We also did not reach the precision wanted for the tossing scenario when objects are picked with a suction cup. First, contact location uncertainties can occur, but this can be mitigated with online inertia identification (including the CoM). The challenging issue is the non-linear deformation of the suction cup. Also, a drop of pressure at release induces non-negligible perturbations. It is crucial to integrate the suction-cup dynamics in the planning as an extension of the robot model, e.g., [START_REF] Cheng | Manipulation with suction cups using external contacts[END_REF], [START_REF] Pham | Critically fast pick-and-place with suction cups[END_REF]. Another issue is the uncertainty in the release moment that can lightly change the object's ballistic trajectory. Since the latter depends on the release pose and velocity, the target can become out of reach. We also noticed that extra care should be given to the end-effector orientation motion at release to prevent undesired interaction between the object and the end-effector after it is released. The experiments also confirmed that for boxing the impacting location and impacting velocities are very sensitive parameters for the success of the tilting.

VIII. CONCLUSION AND FUTURE WORK

We propose a decoupled approach to plan motions for a class of impact tasks based on kinodynamic sampling in jointspace. The new planner is plugged into modules for tossing, grabbing, and boxing, each of which links its target impact objective to via-point configuration parameters (in terms of joint-space position and velocity), which our planner takes into account. The entire task is executed by a standard taskspace QP controller, which is proven to be highly efficient for generating intentional impacts. This work will be further investigated by benchmarking our developed methods together with others (data driven [START_REF] Bombile | Dual-arm control for coordinated fast grabbing and tossing of an object: Proposing a new approach[END_REF], [START_REF] Khurana | Learning to hit: A statistical dynamical system based approach[END_REF] and reference spreading [START_REF] Van Steen | Dual arm impact-aware grasping through time-invariant reference spreading control[END_REF]) in real logistic cases that include robotics with our industrial partner in the frame of the I-AM project. We will also investigate the possibility of casting the entire tasks as a global optimization problem, with more robustness to uncertainties by incorporating suction cup models and additional tasks such as impact assemblies and on-the-fly catching of objects.

Algorithm 2 1 :Fig. 3 :

 213 Fig. 3: Planning framework flowchart.

Fig. 6 :

 6 Fig. 6: Dual grabbing scenario: (left) define impact-grab velocity, (right) execute the planned trajectory.

7

 7

  .

Fig. 7 :

 7 Fig.7: On the left side, the force applied on the box exactly before starting to tilt. On the right one, the box and applied forces during tilting phases S 1 and S 2 .

Fig. 8 :

 8 Fig. 8: Position, velocity, and acceleration tracking for joint 1: desired (dark) and real (shaded), (Top-left). Position and velocity tracking errors (Top-right). Measured Position and Velocity errors to the TS defined Tossing Points (Bottom).

Fig. 9 :

 9 Fig. 9: Same target with impact velocities.

This paper was recommended for publication by Editor Chao-Bo Yan upon evaluation of the Associate Editor and Reviewers' comments. This work was supported by the Research Project I.AM. through the European Union H2020 program (GA 871899) A. Zermane and A.

e-commerce sorting logistics are cardboard boxes; and we used boxes obtained from real deliveries.

Figure 10 shows the execution of the planned trajectory. The desired impact point M is reached with the corresponding velocity to induce the intended planned rotational motion (all obtained from computations in Sec. VI). q (rad/s)


Keywords: A: titanium, A: alloy, B: IR spectroscopy, B: XRD, B: SEM, C: high temperature corrosion

The oxidation kinetics of Ti-6Al-4V samples nitrided in forming gas (N 2 +5%H 2 ) between 520 and 805 • C is investigated using infrared emissivity measurements. Several kinetic exponents have been observed depending on temperature and time. The activation energy corresponding to the main parabolic rutile growth stage is 223 ± 10 kJ/mol. This value is significantly higher than that corresponding to TiN oxidation, which suggests that it is not the rate-limiting step. The protective role of the nitridation is attributed to the internal Ti 2 N layer. The microstructural characterization reveals a reduction in the oxygen-diffusion depth compared to non-nitrided Ti-6Al-4V, and a multi-layered Al 2 O 3 -TiO 2 oxide.

Isothermal oxidation kinetics of nitrided Ti-6Al-4V

studied by infrared emissivity

Introduction

The high strength-to-weight ratio and excellent biocompatibility of titanium and its alloys have been widely exploited in sectors such as aeronautics and biomedicine. In particular, the α + β Ti-6Al-4V alloy remains the world standard for these applications, traditionally accounting for the largest share of the global titanium alloy production [START_REF] Boyer | An overview on the use of titanium in the aerospace industry[END_REF][START_REF] Cui | Titanium alloy production technology, market prospects and industry development[END_REF]. Initially developed for aerospace and aircraft components, it is now a staple of marine applications, power generation, automotive parts, sports, and medicine, with further applications expected to grow even further due to improvements in processing technology [START_REF] Liu | Additive manufacturing of Ti6Al4V alloy: A review[END_REF].

The main disadvantage of titanium alloys is its poor resistance to hightemperature corrosion and wear, which limit its working temperature range to about 400 • C and can also be troublesome in joint implants, where friction may lead to leaching of dangerous ions to the body. A number of treatments have been attempted to improve on these surface properties, with thermal oxidation (TO) being considered one of the simplest and most cost-effective [START_REF] Lin | Surface damage mitigation of titanium and its alloys via thermal oxidation: A brief review[END_REF][START_REF] Maestro | Cyclic thermal oxidation evaluation to improve Ti6Al4V surface in applications as biomaterial[END_REF][START_REF] Kumar | Thermal oxidation of Ti6Al4V alloy: Microstructural and electrochemical characterization[END_REF][START_REF] Guleryuz | Effect of thermal oxidation on corrosion and corrosion-wear behaviour of a Ti-6Al-4V alloy[END_REF].

TO processes require good knowledge of corrosion kinetics and morphology. Moreover, these data are also essential if the working temperatures of this alloy are to be increased. However, literature data about this topic reveal a generally complex oxidation behaviour, with kinetic laws ranging from linear to cubic at different temperatures and time scales, as well as differences in activation energy and morphology [START_REF] Du | Air oxidation behaviour of Ti-6Al-4V alloy between 650 and 850[END_REF][START_REF] Chaze | Influence of aluminium on the oxidation of titanium between 550 and 750 • C[END_REF][START_REF] Frangini | Various aspects of the air oxidation behaviour of a Ti6Al4V alloy at temperatures in the range 600 -700 • C[END_REF][START_REF] Guleryuz | Oxidation of Ti-6Al-4V alloy[END_REF].

Nitriding treatments have been known for a long time to be useful for improving the wear resistance and other tribological properties of titanium alloys [START_REF] Bell | Surface engineering of titanium with nitrogen[END_REF][START_REF] Spies | Second Lakhtin memorial lecture[END_REF]. Among them, gas nitriding has also been suggested as a potentially effective treatment against oxidation [START_REF] Pérez | Influence of nitriding on the oxidation behaviour of titanium alloys at 700 • C[END_REF][START_REF] Lee | Gas nitriding and subsequent oxidation of Ti-6Al-4V alloys[END_REF][START_REF] Lee | Gas nitriding and oxidation of Ti-6Al-4V alloy[END_REF]. However, to our knowledge, only one qualitative study of the temperature-dependent oxidation kinetics of nitrided Ti-6Al-4V has been performed [START_REF] Lee | Oxidation of nitride layers formed on Ti-6Al-4V alloys by gas nitriding[END_REF].

This work contains a study on the oxidation kinetics and morphology of Ti-6Al-4V samples pre-nitrided in N 2 +5%H 2 using infrared emissivity measurements. Conventional oxidation kinetics studies rely mostly on weight-based measurement methods. Alternative techniques capable of determining oxidation kinetics include those based on optical techniques, such as emissivity measurements. This procedure takes advantage of the spontaneous thermal radiation emitted by the sample itself at the high temperatures required for oxidation. The use of alternative techniques (mainly those capable of higher sampling frequencies) is useful for reducing systematic errors and cross-validating results. The method followed in this work has been described previously in the literature [START_REF] Del Campo | Iron oxidation kinetics study by using infrared spectral emissivity measurements below 570 • C[END_REF][START_REF] Del Campo | Kinetics inversion in isothermal oxidation of uncoated WC-based carbides between 450 and 800 • C[END_REF]. Studies on the evolution of the emissivity during oxidation of both titanium [START_REF] Neuer | Insitu measurements of layer thickness during oxidation of titanium[END_REF][START_REF] Zhang | Infrared spectral emissivity property of pure titanium in the 473-1035 K temperature range[END_REF] and Ti-6Al-4V [START_REF] Li | Study of Ti-6Al-4V alloy spectral emissivity characteristics during thermal oxidation process[END_REF] have been made, but the discussion of the oxidation kinetics and morphology has been scarce and qualitative. The overall purpose of this study is to evaluate the effects of nitriding in an N 2 /H 2 atmosphere for the subsequent thermal oxidation of Ti-6Al-4V.

Experimental details

The samples used in this study were cut from the same batch of commercial Ti-6Al-4V (6.59% Al, 4.10 % V, 0.24 % Fe, 0.002 % N, 0.21 % O). The samples were ground with 80 to 320 grit SiC paper and cleaned in an acetone ultrasonic bath. The roughness of the samples was measured with a conventional roughness tester before the treatments. A range of values from 0.25 to 0.48 µm was found for the average roughness of the five samples (R a = 0.38 ± 0.09 µm). This roughness level is similar to others in the literature [START_REF] Lee | Gas nitriding and subsequent oxidation of Ti-6Al-4V alloys[END_REF][START_REF] Lee | Gas nitriding and oxidation of Ti-6Al-4V alloy[END_REF][START_REF] Lee | Oxidation of nitride layers formed on Ti-6Al-4V alloys by gas nitriding[END_REF].

Oxidation kinetics were determined using the infrared emissivity method [START_REF] Del Campo | Iron oxidation kinetics study by using infrared spectral emissivity measurements below 570 • C[END_REF][START_REF] Del Campo | Kinetics inversion in isothermal oxidation of uncoated WC-based carbides between 450 and 800 • C[END_REF]. Time-dependent emissivity measurements were carried out using the high-accuracy infrared emissometer HAIRL [START_REF] Del Campo | New experimental device for infrared spectral directional emissivity measurements in a controlled environment[END_REF], formed by three main modules: a vacuum sample chamber, a FTIR spectrophotometer and a blackbody source for radiance calibration. The vacuum chamber allows to carry out the measurements in a controlled atmosphere. The sample temperature is measured using two bare Type K thermocouples spot-welded on the sample surface out of the area viewed by the detector. The emissivity is measured directly using the blacksur method [START_REF] Pérez-Sáez | Analysis of the accuracy of methods for the direct measurement of emissivity[END_REF] and the modified two-temperature calibration method [START_REF] González-Fernández | Analysis of calibration methods for direct emissivity measurements[END_REF]. All the sources of emissivity uncertainty have been analyzed for this apparatus [START_REF] Del Campo | Combined standard uncertainty in direct emissivity measurements[END_REF]. A spectral resolution of 16 cm -1 was chosen, with wavelength accuracy being ensured by using an FTIR spectrophotometer.

The emissivity results acquired during the oxidation process of the five samples can be employed to determine the oxidation kinetics. As the oxide scale grows, the emissivity spectra show interferential effects due to reflections at the interfaces [START_REF] Del Campo | Iron oxidation kinetics study by using infrared spectral emissivity measurements below 570 • C[END_REF][START_REF] Del Campo | Kinetics inversion in isothermal oxidation of uncoated WC-based carbides between 450 and 800 • C[END_REF]. The oxide scale thickness is related to the maxima and minima position in the spectra according to the following equations:

d = 2m -1 4 λ (1)
d = m 2 λ ( 2
)
where m stands for the order of the interference, λ is the wavelength of the interferential maximum (1) or minimum (2) and d (hereafter simply d) is the apparent oxide film thickness (real thickness multiplied by the real part of the refractive index of the oxide) [START_REF] Del Campo | Kinetics inversion in isothermal oxidation of uncoated WC-based carbides between 450 and 800 • C[END_REF]. Eight samples were used in this study. Three samples were used to test the effects of hydrogen content in the nitridation of this alloy: one acted as a control sample, while the two others were treated at 850 • C for 2 h in 1-bar atmospheres of high-purity N 2 (≤ 3 ppmv H 2 O, ≤ 2 ppmv O 2 , ≤ 0.5 ppmv hydrocarbons) and commercial forming gas (N 2 +5%H 2 ), respectively. Then, the five remaining samples were treated in N 2 +5%H 2 at 850 • C for 2 h and then isothermally oxidized at different temperatures. Measurements of the normal spectral emissivity were used as indicators of the surface state of the sample during the process, as it was observed that it is sensitive to changes in the surface [START_REF] González-Fernández | Infrared normal spectral emissivity of Ti-6Al-4V alloy in the 500-1150 K temperature range[END_REF]. After the treatment was finished, the samples were stabilized at the required temperature for oxidation, and the chamber was opened to laboratory air. The normal spectral emissivity was then measured in a continuous mode during the corrosion process.

Phase identification was carried out by X-ray diffraction (XRD) measurements in the Bragg-Brentano geometry with Cu K α radiation and a step increment of 0.026 • . The surface and cross-section morphologies and distributions of the elements were studied using a JEOL JSM-6400 scanning electron microscope (SEM) equipped for energy-dispersive X-ray spectroscopy (EDX). Images were taken at a working distance of 15 mm with backscattered electrons for compositional contrast and an accelerating voltage of 20 kV. EDX maps correspond to the distributions of the intensities of K α lines for all components except for V, for which the K β lines have been recorded to avoid overlap with Ti lines.

Results and discussion

Nitridation

Fig. 1 shows the XRD diffractograms of the first set of samples used to determine the influence of nitridation in the two atmospheres under consideration. The pattern corresponding to the untreated sample features peaks corresponding to the α and β phases. Significant differences are observed in the case of the nitrided samples, where both TiN and Ti 2 N phases can be observed, with higher integrated intensities in the case of the sample nitrided in forming gas.

This could be expected due to the catalytic effect of hydrogen in nitridation, which allows a deeper penetration of nitrogen in the titanium matrix [START_REF] Spies | Second Lakhtin memorial lecture[END_REF].

A second observation concerns the evolution in the α phase of the nitrided samples. Some of its peaks experience a noticeable shift to lower angles. This corresponds to an increase in the c/a ratio [START_REF] Malinov | Synchrotron X-ray diffraction study of the phase transformations in titanium alloys[END_REF], which has been shown to correlate strongly to the hardness of the material [START_REF] Ashrafizadeh | Structural features and corrosion analysis of thermally oxidized titanium[END_REF][START_REF] Borgioli | Glow-discharge and furnace treatments of Ti-6Al-4V[END_REF][START_REF] Yan | Surface hardening of titanium by thermal oxidation[END_REF]. It is thus shown that the nitriding treatment hardens the subsurface in a similar way for both atmospheres. In the case of the (002) peak at 38 -39 • , a shoulder to the left of the main peak appears for both nitrided samples. A similar feature seems to appear in the following (101) peak at 40 -41 • for the forming-gas-treated sample, whereas the feature in the N 2 -treated diffractogram can be better accounted for by some remaining cubic β phase (with a shortened cell parameter with respect to the untreated sample). The presence of similar shoulders has been noted in the literature [START_REF] Malinov | Synchrotron X-ray diffraction study of the phase transformations in titanium alloys[END_REF][START_REF] Sofinowski | In situ characterization of a high work hardening Ti-6Al-4V prepared by electron beam melting[END_REF]. These minor peaks hint at an α-like phase with similar cell parameters than both the original α and the hardened α, which is sometimes tentatively described as martensite α . Their stronger presence for the (002) peak (38 -39 • ) is perhaps consistent with the larger breadth of this peak for the untreated sample, which is linked to the presence of residual martensite in the original sample [START_REF] Lima | The influence of nitrogen to hydrogen ratio and temperature on thickness and phase composition in plasma nitrided Ti-6Al-4V[END_REF].

Oxidation kinetics

Five samples have been oxidized at temperatures from 520 • C to 805 • C. Fig. 2 shows, as an example, the interferential emissivity oscillations as a function of oxidation time for various wavelengths at 600 and 740 • C. In agreement with Eq. ( 1), in both cases the time needed to reach the first-order interferential maximum increases with wavelength. It is worth mentioning that the time required to observe the first maximum at 600 • C is almost two orders of magnitude higher than the one for 740 • C. These interferential patterns stop when the oxide layer is thick enough to become opaque at such wavelengths. This highlights the advantage of a wide spectral range when using optical techniques for studying oxidation kinetics, with the mid-infrared range being the one that remains semi-transparent for the longest time [START_REF] González De Arrieta | Evolution of the infrared emissivity of Ni during thermal oxidation until oxide layer opacity[END_REF].

The time dependence of the apparent oxide thickness can be obtained directly from the interferential maxima and minima by applying Eqs. ( 1) and (2). Following Ref. [START_REF] Du | Air oxidation behaviour of Ti-6Al-4V alloy between 650 and 850[END_REF], log-log plots of the growth of the oxide layer thickness with time at each temperature are shown in Fig. 3. The kinetic exponents have been calculated using the standard empirical equation:

d n = kt ( 3
)
where d is the apparent oxide thickness, n the kinetic exponent, k the apparent rate constant and t the oxidation time.

In the case of the samples oxidized at temperatures below 700 • C (Figs. 3ac), the oxidation process takes place with a single well-defined exponent, except for the early stages of the oxidation at 520 • C, where a transitory process of linear-parabolic nature is observed. This is consistent with a surface mode of oxidation transforming into a diffusion process once a sufficiently thick layer is formed. It is important to mention that this method is capable of detecting apparent oxide layer thicknesses of at least one-quarter of the lowest measurable wavelength. Therefore, it cannot probe the earliest stages of the oxidation below such thicknesses, and thus it cannot discard a similar transitory stage for the other oxidation temperatures (such as the one hinted at in Fig. 3c). These results show a clear temperature dependence of the kinetic exponent, decreasing clearly from near-cubic (2.58) to near-parabolic (1.86) to linear-parabolic (1.69). Despite this evolution, the overall kinetics can still be described within the same framework (corresponding to the diffusive growth of a rutile layer), similarly to the results reported in Ref. [START_REF] Frangini | Various aspects of the air oxidation behaviour of a Ti6Al4V alloy at temperatures in the range 600 -700 • C[END_REF].

Figure 3 also shows that, in the case of temperatures above 700 • C, three processes with similar structure can be observed. The oxidation follows firstly a clear parabolic behaviour which is interrupted by a temporary diversion into cubic kinetics, before proceeding on to the final stage, which is either parabolic again (740 • C, Fig. 3d) or linear (805 • C, Fig. 3e). This is consistent with the observed spallation in the case of the sample oxidized at 805 • C, and also with literature data, which sets 800 • C as the limit of parabolic kinetics [START_REF] Du | Air oxidation behaviour of Ti-6Al-4V alloy between 650 and 850[END_REF]. It should be noted that the oxidation process shown in this article at 805 • C does feature higher-than-linear kinetics for the first three hours of oxidation (n = 1.93), which suggests a regular diffusion-limited oxidation process. Other works have also reported such type of kinetics, with spallation occurring only above a certain time threshold [START_REF] Du | Air oxidation behaviour of Ti-6Al-4V alloy between 650 and 850[END_REF]. This indicates that the linear kinetics around 800 • C are not due to an intrinsically defective oxide layer but to a susceptibility to crack and the inability of the oxide layer to recover at high temperatures. Overall, the observed kinetics for the oxide layer growth is complex but consistent with the better-established cases of non-nitrided Ti-6Al-4V and Ti.

Activation energy

The main parabolic or quasi-parabolic (in the case of 520 • C and 665 • C) stage of each of the five temperatures in Fig. 3, which are assumed to correspond to the growth of the main rutile layer, have been fitted to the following formula [START_REF] Pieraggi | Calculations of parabolic reaction rate constants[END_REF]:

(d -d 0 ) 2 = k p (t -t 0 ) (4)
where k p corresponds to the parabolic constant, d 0 is the initial thickness, and t 0 is the duration of the faster transient oxidation stage. An activation energy of 223 ± 10 kJ/mol has been found, as shown in Fig. 4. This energy is similar to that obtained for pure titanium [START_REF] Chaze | Influence of aluminium on the oxidation of titanium between 550 and 750 • C[END_REF] and higher than some reported values for Ti-6Al-4V [START_REF] Frangini | Various aspects of the air oxidation behaviour of a Ti6Al4V alloy at temperatures in the range 600 -700 • C[END_REF], but much lower than those reported for the corresponding parabolic kinetics in Ti-6Al-4V and Ti-Al alloys (∼ 270 kJ/mol) [START_REF] Du | Air oxidation behaviour of Ti-6Al-4V alloy between 650 and 850[END_REF][START_REF] Chaze | Influence of aluminium on the oxidation of titanium between 550 and 750 • C[END_REF][START_REF] Guleryuz | Oxidation of Ti-6Al-4V alloy[END_REF]. It should be mentioned that the activation energies obtained through weight-gain measurements account for the total oxygen uptake (i.e., oxide layer growth plus oxygen diffusion), whereas those based in thickness measurements account only for the growth of the oxide layer. Out of the aforementioned works, only Ref. [START_REF] Frangini | Various aspects of the air oxidation behaviour of a Ti6Al4V alloy at temperatures in the range 600 -700 • C[END_REF] reports a separate activation energy for oxide layer growth in air for Ti-6Al-4V (147 kJ/mol), which is significantly lower than the value reported in this work. The reported value is therefore much closer to the self-diffusion activation energy of rutile (251 kJ/mol) [START_REF] Young | High temperature oxidation and corrosion of metals[END_REF]. Finally, the observed activation energy is also significantly higher than the values reported for TiN (110-136 kJ/mol), which suggests that the oxidation of this phase is not the rate-limiting step in this case [START_REF] Chen | Oxidation behavior of titanium nitride films[END_REF]. On the contrary, Ti 2 N is known to be a better barrier against diffusion [START_REF] Abdallah | STEM-EELS identification of TiO X N Y , TiN, Ti 2 N and O, N dissolution in the Ti2642S alloy oxidized in synthetic air at 650 • C[END_REF][START_REF] Pohrelyuk | Effect of oxidation of nitride coatings on corrosion properties of Ti-6Al-4V alloy in 0.9% NaCl at 40 • C[END_REF], although no activation energy has been found for the oxidation of this compound in air. Interestingly, it has been established that oxide growth on top of nitride phases proceeds by outward titanium diffusion [START_REF] Pérez | Influence of nitriding on the oxidation behaviour of titanium alloys at 700 • C[END_REF], contrary to the inward oxygen diffusion that is mostly responsible for the oxidation of pure titanium [START_REF] Du | Air oxidation behaviour of Ti-6Al-4V alloy between 650 and 850[END_REF]. This may represent a fundamental difference in the oxidation process, which remains throughout the rutile growth stage. This may be due to the internal Ti 2 N layer controlling the nature of the diffusing species.

It has been claimed that nitridation at high temperatures leads to the same oxidation rates for pure titanium and Ti-6Al-4V [START_REF] Pérez | Influence of nitriding on the oxidation behaviour of titanium alloys at 700 • C[END_REF]. It has also been stated that rutile formed in the presence of nitrogen features a larger unit cell volume than pure TiO 2 , due to the substitution of O ions by trivalent N ions. The presence of defects alters the diffusion paths available for O, and could, therefore, affect the activation energy [START_REF] Borbunov | The role of nitrogen in the oxidation of titanium in air at elevated temperatures[END_REF]. Thus, a tentative hypothesis may involve a negation of the beneficial effects of the Al content when the concentration of N ions in rutile is so large. Moreover, the concentration of Al in the surface is reduced due to its low solubility in titanium nitrides and the observed segregation in nitrided Ti-6Al-4V [START_REF] Lee | Gas nitriding and subsequent oxidation of Ti-6Al-4V alloys[END_REF]. However, the absence of more data on the specific activation energy associated with the growth of the oxide layer itself makes comparisons to regular Ti-6Al-4V difficult. To our knowledge, no activation energy data on nitrided Ti-6Al-4V is available in the literature, as most literature on the topic has focused on the microstructure of the nitrided alloy, and not on the kinetics of the process. Even when kinetics have been studied [START_REF] Pérez | Influence of nitriding on the oxidation behaviour of titanium alloys at 700 • C[END_REF][START_REF] Lee | Oxidation of nitride layers formed on Ti-6Al-4V alloys by gas nitriding[END_REF], there have been no reports on activation energy values.

XRD characterization of the oxidized samples

Fig. 5 shows the XRD patterns corresponding to samples oxidized at low, medium and high temperatures. Rutile and alumina are the only oxidation products seen, while a still significant Ti 2 N presence can also be observed for the sample oxidized at 520 • C. This residual presence, together with the quasi-cubic kinetics, suggests that the formation of nitrides is effective against corrosion at relatively low temperatures. It should be noted that no anatase TiO 2 has been found for any sample, in contradiction with a number of studies on non-nitrided Ti-6Al-4V [START_REF] Kumar | Thermal oxidation of Ti6Al4V alloy: Microstructural and electrochemical characterization[END_REF][START_REF] Guleryuz | Effect of thermal oxidation on corrosion and corrosion-wear behaviour of a Ti-6Al-4V alloy[END_REF][START_REF] Guleryuz | Oxidation of Ti-6Al-4V alloy[END_REF][START_REF] Li | Study of Ti-6Al-4V alloy spectral emissivity characteristics during thermal oxidation process[END_REF][START_REF] Ashrafizadeh | Structural features and corrosion analysis of thermally oxidized titanium[END_REF]. As anatase is less protective than rutile against certain modes of corrosion (e.g., by reducing acids) [START_REF] Guleryuz | Effect of thermal oxidation on corrosion and corrosion-wear behaviour of a Ti-6Al-4V alloy[END_REF], its absence after nitriding could be positive for the development of this method of protection. On the other hand, improved osteogenic activities have been suggested for mixed anataserutile layers in pure Ti [START_REF] Wang | Surface thermal oxidation on titanium implants to enhance osteogenic activity and in vivo osseointegration[END_REF]. Oxidation of nitrided samples in the literature was performed at temperatures higher than those where anatase is expected to form, so no comparison is possible. Nevertheless, it has been stated that the amount of anatase formed at 500 • C for short oxidation times decreases with time at the same time that the amount of Ti(O) phase grows [START_REF] Kumar | Thermal oxidation of Ti6Al4V alloy: Microstructural and electrochemical characterization[END_REF].

In order to focus on the evolution of the α phase (in connection with Section 3.1), X-ray diffractograms of the samples oxidized at the lowest temperatures are shown in Fig. 6. Interestingly, peaks corresponding to Ti 2 N (but not TiN) are observed in the diffractogram corresponding to the sample oxidized at 520 • C, which indicates that this phase remains during the rutile-forming stages. Besides, a decrease in the amount of α phase is noted, together with a surge in a new phase with the same space group and similar lattice constants. This αlike phase, which was already observed in Section 3.1, has not been conclusively identified. Similar shoulders have been reported in the literature, mostly linked to the hcp martensite α' phase [START_REF] Malinov | Synchrotron X-ray diffraction study of the phase transformations in titanium alloys[END_REF][START_REF] Sofinowski | In situ characterization of a high work hardening Ti-6Al-4V prepared by electron beam melting[END_REF]. Precipitation of α 2 -Ti 3 Al in the subsurface of oxidized Ti-6Al-4V has also been described in [START_REF] Bagot | The effect of oxidation on the subsurface microstructure of a Ti-6Al-4V alloy[END_REF]. Other complex subsurface microstructures have been reported, for example, for oxidized Ti-6Al-4V-1B, in which the presence of an additional element (boron) stimulates the precipitation of new phases [START_REF] Brice | Oxidation behavior and microstructural decomposition of Ti-6Al-4V and Ti-6Al-4V-1B sheet[END_REF]. The information shown in this angular range is supported by extra peaks in other ranges, but no ultimate conclusion can be inferred from the data to suggest a definitive indexation. In any case, the structure of the α phase is not critical for the interpretation of this work.

Similar to the analysis in Section 3.1, the angular shift of the α peaks can be correlated with an increase in the c/a ratio. A similar shoulder to the one found for the nitrided samples can be seen for the sample oxidized at the lowest temperature (520 • C), whereas the shape of the diagram for the sample oxidized at 600 • C hints at a double peak feature for the (002) peak at 38 • . Finally, multiple α-like phases are no longer observed for the sample oxidized at 665 • C. Interestingly, even the (100) peak at 35 -36 • shifts to lower angles in this case. This corresponds to a change in all the lattice parameters, not only the c/a ratio, as this peak depends only on the a parameter. It can be seen that oxidation at low temperatures does not induce significant changes to the subsurface of the alloy with respect to nitridation, but stronger effects are observed at higher temperatures.

SEM characterization of oxidized samples

Surface SEM images for samples oxidized at low, medium and high temperatures are shown in Fig. 7. The sample oxidized at 520 • C shows the texture of the underlying substrate, with very little oxidation. Small nucleation sites with sizes much lower than 1 µm can be observed. As temperature rises, the formation of granular oxide particles can be observed, which cover the surface completely, hiding the substrate morphology. For the oxidation at 665 • C (Fig. 7b), the particle size distribution is relatively homogeneous at 0.5-1.5 µm, whereas for the highest temperature (Fig. 7c) a more irregular morphology is observed, with significant protuberances and more heterogeneous particle sizes (ranging from 0.25 to 2-3 µm).

The surface morphology of all samples is very similar to that in Ref. [START_REF] Pérez | Influence of nitriding on the oxidation behaviour of titanium alloys at 700 • C[END_REF], where it was shown that nitridation at high temperatures led to smoother oxidized surfaces with finer and better-crystallized grains than in the non-nitrided case, independent of the underlying microstructure. In particular, the formation of coarse nodules is avoided by nitridation. This may be of interest for biomedical applications, as sharper and more homogeneous surface morphologies show improved wettability [START_REF] Wang | Surface thermal oxidation on titanium implants to enhance osteogenic activity and in vivo osseointegration[END_REF].

A cross-sectional SEM image of the sample oxidized at 805 • C is shown in Fig. 8. As seen in other references for nitrided Ti-6Al-4V [START_REF] Pérez | Influence of nitriding on the oxidation behaviour of titanium alloys at 700 • C[END_REF], there is a relatively large amount of β phase in the subsurface of the alloy, contrary to the more established case of regular Ti-6Al-4V [START_REF] Du | Air oxidation behaviour of Ti-6Al-4V alloy between 650 and 850[END_REF][START_REF] Chaze | Influence of aluminium on the oxidation of titanium between 550 and 750 • C[END_REF][START_REF] Frangini | Various aspects of the air oxidation behaviour of a Ti6Al4V alloy at temperatures in the range 600 -700 • C[END_REF][START_REF] Guleryuz | Oxidation of Ti-6Al-4V alloy[END_REF]. Whereas the concentration of subsurface β phase is lower than in the bulk of the material, it is still present, and so no α case can be said to be formed for this alloy. Despite the significant evolution of the α phase evidenced by XRD, there is no formation of a β-less case several µm deep. This effect of nitridation has been corroborated in other works [START_REF] Pérez | Influence of nitriding on the oxidation behaviour of titanium alloys at 700 • C[END_REF]. Nevertheless, it should be noted that, also because of nitridation, the hardening effect shown in Fig. 6 is taking place in a much thinner diffusion zone, when compared with the standard results of oxidation in air. This could be advantageous for improving the mechanical properties of the alloy. Hardness tests have not been attempted due to the small size of the diffusion layer and the significant presence of β grains near the surface, which makes the hardness inhomogeneous at the µm range [START_REF] Pérez | Influence of nitriding on the oxidation behaviour of titanium alloys at 700 • C[END_REF].

In order to discern the origin of the cubic oxidation regimes at such temperatures, a compositional EDX map of the oxide layer is shown in Fig. 9 (a different location is shown in the Supplementary Material). It can be easily seen that the oxidation product is not a uniform layer, but consists of a top alumina layer, and two thick rutile sublayers separated by a much thinner alumina layer located below a lateral crack. Such morphologies can be observed in some reported micrographs of oxidized Ti-6Al-4V, although not as the focus of the study [START_REF] Guleryuz | Oxidation of Ti-6Al-4V alloy[END_REF]. Based on the observation of a periodic stacking of rutile and alumina layers, a model was formulated for the oxidation of Ti-6Al-4V [START_REF] Du | Air oxidation behaviour of Ti-6Al-4V alloy between 650 and 850[END_REF], in which these bilayers grow again after each crack of the previous layer. The formation of alumina between two rutile growth stages is consistent with the cubic rates (or higher) found for the samples oxidized at the highest temperatures (Figs. 3d-e). This is also in accordance with other experimental results, which indicate that the growth of Al 2 O 3 in alloys generally follows such kind of exponents [START_REF] Tallman | A critical review of the oxidation of Ti 2 AlC, Ti 3 AlC 2 and Cr 2 AlC in air[END_REF].

Conclusions

The effects of nitriding in forming gas in the oxidation kinetics and morphology of Ti-6Al-4V have been studied. It is shown that the hydrogen content allows for more effective formation of nitrides in this alloy. The presence of such layers is advantageous, leading to finer surface morphologies, thinner diffusion layers, absence of anatase and mostly protective oxide layers. The kinetics of the oxidation process has been discussed in terms of the mechanisms behind them. The samples oxidized at temperatures below 700 • C feature a quasi-parabolic behaviour, with a rate exponent which steadily decreases with temperature. In the case of the two highest temperatures, a more complex behaviour has been observed, due to the cracking of the rutile layer and the formation of an intermediate alumina layer. An activation energy of 223 ± 10 kJ/mol has been found for the rutile forming stages. Comparison of this value to literature data for non-nitrided Ti-6Al-4V is difficult, as reported values usually correspond to the total oxygen intake, instead of the oxide layer growth.

2 Figure 1 :

 21 Figure 1: XRD diffractograms of the untreated sample and those treated in N 2 +5%H 2 and N 2 .

Figure 2 :

 2 Figure 2: Temporal evolution of the normal spectral emissivity of Ti-6Al-4V at (a) 600 • C and (b) 740 • C.

Figure 3 :

 3 Figure 3: Growth of the oxide layer apparent thickness with time at different oxidation temperatures. Kinetic exponents obtained by fitting to Eq. (3) are shown for each distinct oxidation stage.

Figure 4 : 3 Figure 5 :

 435 Figure 4: Arrhenius plot of the rutile-forming parabolic oxidation stages obtained by fitting the quasi-parabolic regions to Eq. 4. Activation energy: 223 ± 10 kJ/mol.

Figure 6 :

 6 Figure 6: Close-up X-ray diffractograms of the three samples oxidized at lower temperatures.

  (a) 520 • C. (b) 665 • C. (c) 805 • C.

Figure 7 :

 7 Figure 7: Surface SEM image of the oxide layer morphology at different temperatures.

Figure 8 :

 8 Figure 8: Large-area cross-sectional SEM image of the sample oxidized at 805 • C. The presence of β grains (light colour) near the surface is observed.

  (a) Cross-sectional SEM image of the oxide layer. (b) EDX spectrum of Ti. (c) EDX spectrum of Al. (d) EDX spectrum of O. (e) EDX spectrum of V

Figure 9 :

 9 Figure 9: Cross-sectional SEM image and compositional EDX maps of the the substrate-oxide layer interface.

Acknowledgements

This work has been carried out with the financial support of the SAIOTEK program (Project number S-PC08UN07) of the Basque Government and the Universidad-Empresa program (Project number UE06/01) of the University of the Basque Country UPV/EHU, in collaboration with Industria de Turbo Propulsores S.A. I. González de Arrieta would like to acknowledge the Basque Government for its support through a PhD fellowship (PRE-2018-1-0045 and PRE-2019-2-0191). The authors of this paper recognize and thank the leadership in infrared emission spectroscopy of the recently deceased Prof. Raúl B. Pérez-Sáez, whose last work is this contribution.

Data statement

Data will be made available on request.

Supplementary material


Nabila Mansouri
email: nabila.elmansouri@gmail.com

Yousra Ben Jemaa
email: yousra.benjemaa@enis.rnu.tn

Cina Motamed
email: motamed@lisic.univ-littoral.fr

Antonio Pinti
email: antonio.pinti@univ-valenciennes.fr

Eric Watelain
email: eric.watelain@univ-valenciennes.fr

A NEW STRATEGY BASED ON SPATIOGRAM SIMILARITY ASSOCIATION FOR MULTI-PEDESTRIAN TRACKING

Keywords: Multi-pedestrian tracking, Spatiogram, Coupling detection, Tree structure, Interpolation, Human detector

Multiple pedestrian tracking is an active and challenging research topic that many different approaches have addressed it. Since Human stick changes over time and person usually moving in random way, the identity association remains a hard task. In this paper we propose a new method for coupling detections over all the frames of the video sequence in order to make a performant tracker. The Pedestrian detection is ensured using the Dalal and Triggs's human detector. In order to overcome the problem of missing detections caused by occlusion, we propose to use an interpolation process based on average speed. Then, all the previously detections are organized over a tree structure, where each frame represents a tree level. All detections in level 'i' are linked to the next level by an arc characterized by a cost representing the spatiogram similarity between these 2 detections. After trajectories refinement is done based on Euclidean distance to palliate the false detection association. An experimental study conducted in 2 datasets (CAVIAR and our database) proves the good performance of our proposed method in term of tracker precision and tracker accuracy.

I. INTRODUCTION

Pedestrian tracking is an active research area that can be applied in different fields such as in video surveillance, people behavior analysis, road crossing security. Tracking pedestrian remains a challenging problem in computer vision works given the highly articulated human body. It consists of 2 main stages: 1) pedestrian detection and 2) association of the detected pedestrian windows [START_REF] Yilmaz | Object Tracking: A Survey[END_REF].

1) Pedestrian detection consists on differentiating people from others objects in the image. In this context several works have been proposed in the literature. They can be divided into 3 categories: i) Appearance based approaches, ii) Motion based approaches, iii) Shape based approaches. Works done in [START_REF] Dai | Layered representation for pedestrian detection and tracking in infrared imagery[END_REF], [START_REF] Munder | Pedestrian detection and tracking using a mixture of view-based shape-texture models[END_REF] detect pedestrian by finding the occurrence of skin color [START_REF] Dai | Layered representation for pedestrian detection and tracking in infrared imagery[END_REF] or texture [START_REF] Munder | Pedestrian detection and tracking using a mixture of view-based shape-texture models[END_REF]. Although they are not complex and intuitive, these approaches are very sensitive to occlusion, shadow and lighting variations. Whereas, in [START_REF] Sidenbladh | Detecting human motion with support vector machines[END_REF], [START_REF] Shashua | Pedestrian detection for driving assistance systems : Single-frame classification and system level performance[END_REF] authors use a global or local periodicity of pedestrian motion to discriminate people. Since pedestrian motion is very sensitive to point of view variations and needs a long video sequence in order to detect the periodicity of the displacement, people detection by motion gives unsatisfactory results. The shape approaches are based on a morphological model such as, cylindrical, elliptical, stick models. A probability of similarity is calculated between the predefined model (learning model) and the query one [START_REF] Eng | A bayesian framework for robust human detection and occlusion handling using human shape model[END_REF]. The major limit of this method is the point of view variation between the learning model and the query one that gives a false response. Recently Dalal and Triggs [START_REF] Dalal | Histograms of Oriented Gradients for Human Detection[END_REF] propose a more satisfactory technique to people detection by using the Histogram of Oriented Gradient (HOG) descriptor to represent the target object and a linear Support Vector Machine (SVM) as a classifier. Results given by this approach are insensitive towards color and lighting. But they are sensitive to occlusion, so a pedestrian can be dismissed in some frames. 2) Pedestrian tracking consists on estimating the trajectory of the target previously detected. The motion model like a kalman filter [START_REF] Bertozzi | Pedestrian Localization and Tracking System with Kalman Filtering[END_REF] and particle filter [START_REF] Tung | Human Motion Tracking using a Color-Based Particle Filter Driven by Optical Flow[END_REF] estimates the current state based on the past state recursively. In fact particles filter overcomes the linear motion limitation of kalman filter. But in multi-objects tracking the number of required particles grows considerably and become hardly to do with in practical situations. Also motion models can't handle the trajectory overlap and the identity switch because they haven't any information to discriminate between peoples. For this reason some works propose to merge between the motion and appearance models. The more popular appearance model used in the literature is the color histogram. In fact the approach proposed in [START_REF] Jiang | Multiple Pedestrian Tracking using Colour and Motion Models[END_REF] use kalman filter to predict pedestrian position and 4D color histogram, presented in 'Lab' color space, to ensure identity coupling. Since histogram represents only color distribution, two pedestrians wearing a similar color clothes have an important likelihood to be confused. In [START_REF] P Ã C Rez | Colorbased probabilistic tracking[END_REF], authors propose to combine color histograms defined in 'HSV' space and a particle filter. Given that 'HSV' color space is significantly power compared to others ones the tracker performances are ameliorates. The downside of this method is the important computational cost of the particle filter. The solution adopted by recent works [START_REF] Jiang | Multiple Pedestrian Tracking using Colour and Motion Models[END_REF], [START_REF] Andriyenko | Discrete-Continous Optimization for multi-Target Tracking[END_REF], [START_REF] Kang | Object Reacquisition using Invariant Appearance Model[END_REF] is to associate the previews detections. In fact in [START_REF] Andriyenko | Discrete-Continous Optimization for multi-Target Tracking[END_REF] authors provide a continuous trajectory model based on cubic B-splines. For this purpose, they need to know the start and the end positions. Also to associate inter-frames detections, this method uses only the distance between consecutive positions that remains a considerable shortcoming especially with the trajectories' overlap and (the direction changes). However in [START_REF] Kang | Object Reacquisition using Invariant Appearance Model[END_REF], the authors base their method on the foundations of Bayesian estimation theory. They have only the visibility on the single current object and not on other objects.

In this paper, we search to overcome some challenging problems of multi-pedestrian tracking such as: identity switch, and re-identification after occlusion and pedestrians' trajectories refinement. Our proposed solution is based on coupling detections and data association. Indeed, our main contribution is to re-address the coupling detection process in order to make a multi-pedestrian tracker. So we introduce the tree structure representation of frame by frame detections and exploit the spatiogram descriptor, powerful to ensure a good identity association.

This paper is organized as follows. Section 2 detailed our proposed method based on tracking pedestrian by coupling the detections that maximize the spatiogram similarity, passing by short introduction of Dalel and Trigges' person detection method. Section 3 explains the test protocol and analyzes the performance of our tree research method by presenting some results. Finally, conclusions and future works are discussed in Section 4.

II. COUPLING DETECTIONS BY MAXIMIZING THE SPATIOGRAM SIMILARITY

A. People detection

Our pedestrian tracking process uses Dalal and Triggs's [START_REF] Dalal | Histograms of Oriented Gradients for Human Detection[END_REF] method to detect humans in each video frame. This approach consists of several steps. Firstly, color normalization is done in order to reduce the influence of shadows and illumination changes. The intensity gradient at each pixel is then computed for each color channel. The next step is to acquire a histogram of the gradient orientations at each spatial region or (cell). A range of adjacent cells are grouped into larger blocks so as to normalize edge contrast and illumination. The last stage is to concatenate all of the normalized block vectors over each detection window to form a final window descriptor. Then a linear SVM classifier is applied to successfully classify detection windows into 'non-pedestrian' and 'pedestrian'.

Dalal and Trigss's approach can't detect all pedestrian in each frame because of partial or total occlusion. Indeed we present in the next our contribution in people detection and our coupling process.

B. Coupling detections

1) Spatiogram similarity calculation: Since spatiogram descriptor has good performances according to a comparative study between color, texture and shape descriptors for multipedestrians identification done in [START_REF] Derbel | Comparative study between Color Texture and Shape descriptors for multi-camera pedestrians identification[END_REF], we assume that the use of the spatiogram as an appearance model to check the similarity between targets can improve tracking pedestrian performances. The spatiogram represents an extension of color histogram incorporating spatial information. In this work, we use the Cornaire [START_REF] Conaire | An improved spatiogram similarity measure for robust object localization[END_REF] spatiogram known for its great representational power. First, a quantization phase is performed on the original image to reduce color variation. Second, a color histogram (denoted η b ) was computed, that corresponds to assigning a number of pixels to each color intensity. Then the algorithm calculates for each intensity existing in the quantized image its spatial distribution represented by the mean µ and the covariance σ of the spatial position of all pixels that fall into each histogram bin as explained in equation ( 1)

η b = N i=1 δ ib (1) µb = 1 N j=1 δ jb N i=1 X T i δ ib δ ib = 1 N j=1 δ jb N i=1 (X i -µ b )(X i -µ b ) T δ ib
Where δ ib = 1 if the pixel i falls into bin b and 0 otherwise; C is a normalizing constant ensuring that the sum of all bins is equal to 1 and X i = (x i , y i ) T is the spatial position of each pixel.

To compare two persons spatiograms, S(η, µ, σ) and S (η , µ , σ ) each containing B bins, the following similarity measure presented in equation ( 2) is used.

similarty = B b=1 ψ b η b η b (2)
where

ψ b = ηexp( -1 2 (µ b -µ b )σ -1 b (µ b -µ b )) σ -1 b = σ -1 b + (σ ) -1 b
with σ b represents the spatial similarity measure and η is the Gaussian normalization term.

2) Proposed process for coupling detection: People detection is done frame by frame along the video sequence. Each detection in each frame represents a level in the tree. Every detection in frame i is initially linked to all the detections in frame i + 1 by an oriented arc. The arc cost represents the spatiograms similarity between these 2 linked pedestrians (detections). Spatiograms similarity is calculated as shown previously according to equation [START_REF] Dai | Layered representation for pedestrian detection and tracking in infrared imagery[END_REF]. The pedestrian trajectory is the longest path that connects similar detections. Association of the correct detection to the correct trajectory must ensure the constraint that, in every level the same detection can't be associated to more than one trajectory in the same time.

Given that a frame by frame HOG detection can lost some pedestrians in some frame because of occlusion, our main contribution in people detection in video sequence is to ensure a continuous detection and overcome occlusion problem. In fact, we seek the lost pedestrian position by interpolation process based on average speed. Based on study made in [START_REF] Ismail | Automated Collection of Pedestrian Data Using Computer Vision Techniques[END_REF], walking speed ranging is generally between 3.0 to 4.95 feetper-second. So a pedestrian make approximately an average displacement of 0.05m every video frame. Consequently, we model interpolation process as a translation, with vector -→ u , of the last detected pedestrian windows as shown in figure 1. In fact Bounding box including pedestrian is characterized by a quadruple (x, y, w, h) where (x, y) is the begin point and w and h are respectively the width and height of the bounding box. The translation consists of translating the point P (x, y) to P i (x i , y i ) in the same direction and using a vector -→ u equal to the average distance calculated referred to the average speed by keeping the same dimensions (w and h). Since the false detection will have neither previous nor next similar one. All the coupling detection process is shown in figure 2.

III. PERFORMANCE ANALYSIS

In this section, we present the used datasets and the test protocol. The performance of the proposed system is then presented and analyzed.

A. Used datasets

To evaluate our system performances, we conduct an experimental study on two datasets. The 'CAVIAR' video database [17] having 25 frames/second and 384 × 288 pixels per frame. It integrates 3 peoples walking together along the corridor and sometimes crossing with others with respect to the camera field of view.

The second dataset is ours. We have filmed 2 scenarios. The first contains 2 pedestrians crossing a street on a crosswalk in the same direction and the second contains 3 pedestrians with 2 in the same direction and the third is in opposite direction. These video sequences have 25 frames/second and 720 × 576 pixels per frame. Our database contains voluntarily difficult occlusion situations where people were occulted by vehicles or by other people. Also, it is captured in outdoor context with important lighting changes.

B. Test protocol and results

The test protocol aims to quantify the right coupling detection during the pedestrian displacement. For that, we will introduce some detection coupling in different time of the video sequence. In fact, each pedestrian will be denoted by a unique number (1, 2, 3, ...) during the sequence. The perfect result is found when the pedestrian's associated number remains the same during all the tracking phase. Figure 3 (a) shows the tracking results on CAVIAR database whereas figure 3 (b,c) illustrates the results on our database.

Since the spatiogram performance in person's identification isn't 100% so some false association will happened. In fact in figure 4, blue (and red)points represent respectively the trajectories of the first pedestrian (denoted p1) and the second pedestrian (denoted p2). We remark some blue points in the pedestrian 2 's trajectory, this is caused by some confused identity association. In order to overcome this shortcoming, we conduct a verification process based on natural constraint of human displacement. In fact, a pedestrian walking in the normal condition, can't suddenly change his spatial position.

So the detection point location will be associated to the nearest trajectory if the Euclidean distance between each successive detections associated with the same trajectory bypass the threshold distance as shown in figure 4. The threshold distance is defined as the mean distance between all of the trajectory points. This trajectory filtering gives satisfactory results and corrects the association mistakes. A quantitative evaluation presented above illustrates also the good performance of our trajectories refinement especially in the 'MOTP' value before and after fitting.

C. Quantitative experimental study

An ideal multiple object trackers, should at all points in time, find the correct number of objects present in the video and estimate the position of each object precisely. It should also keep consistent track of each object over time: Each object should be assigned a unique track ID which stays constant throughout the sequence [START_REF] Bernardin | Multiple Object Tracking Performance Metric and Evaluation in a Smart Room Environment[END_REF]. So for a quantitative evaluation we use a widely used metrics proposed by [START_REF] Bernardin | Multiple Object Tracking Performance Metric and Evaluation in a Smart Room Environment[END_REF] such as:

1) The Multiple Object Tracking Precision (MOTP): indicates the coupling precision. MOTP is calculated as explained in equation ( 3)

M OT P = i,t d it t c t (3)
Where d it is the distance between ground-truth trajectories and automatically generated trajectories' positions for each frame t and ct is the total number of matches made in this same frame t . 2) Multi-Object Tracking Accuracy (MOTA): is the average error. It includes the false positives (f pt ), missed targets (m t ) and identity mismatches (mme t ). MOTA is calculated using equation [START_REF] Sidenbladh | Detecting human motion with support vector machines[END_REF].

M OT A = 1 - i,t (f pt , m t , mme t ) t g t (4)
We add also other metric to evaluate the performance of detection process with interpolation technique. The Multiple Object Detection Accuracy (MODA) that checks for missed targets and false positives in each frame t as explained in equation [START_REF] Shashua | Pedestrian detection for driving assistance systems : Single-frame classification and system level performance[END_REF].

M ODA = 1 - i,t (f pt , m t ) t g t (5)
All results are presented in table 1. ing the trajectory refinement. In fact, we ameliorate the trajectory estimation and correct some false associated detections to some trajectories. Our solution is based on the calculus of the Euclidean distance between pedestrian's position in level 'i' and her position in level 'i+1'. If this distance, exceeds the threshold distance the detection will be associated to the nearest trajectory.

A trajectories refinement step ameliorates considerably our tracker precision. In fact the MOTP values are 29.89cm, 10cm respectively before and after trajectories refinement. A decrease of 20cm of the average distance separating the ground truth and estimated trajectories proves the important and the good impact factor of this step.

IV. CONCLUSION AND FUTUR WORKS

In this paper, we have introduced a new approach for multi-pedestrians tracking based on coupling detections. We maintain the good person identity during video sequences by increasing the spatiogram similarity between associated detections. Thanks to our approach, the problem of lost detections caused by occlusion was solved. To ensure this objective, we use an interpolation process based on pedestrian average speed. It is proven that interpolation can proportionally overcomes occlusion impact and gives a good MODA values. Also the detection organization inspired from the tree theory represents a good strategy to perform the global association of pedestrian's identity. An experimental study conducted in 2 datasets (CAVIAR and our database) illustrates the good performances obtained by the proposed method which can reduce the identity confusion problem in multi-pedestrian tracking and gives enhanced fitted trajectories. Future work will focus on developing a robust behavioral analysis module in order to integrate it in the proposed system which use only appearance descriptor (spatiogram). The goal is to improve the global people tracking performances especially over the pedestrian identities estimation in the case of visually similar pedestrians.

Fig. 1 .

 1 Fig. 1. Finding lost detection by interpolation.

Fig. 2 .

 2 Fig. 2. The proposed pedestrian tracking by coupling detections process.

Fig. 3 .

 3 Fig. 3. Tracking results: (a) CAVIAR database, (b) Our dataset: scenario 1, (c) Our dataset: scenario 2.

Fig. 4 .

 4 Fig. 4. Trajectories refinement.

Table 1 .

 1 Quantitative results of the proposed tracking system in term of MOTP, MODA and MOTA.

		MOTP	MODA	MOTA
		Before	After
	CAVIAR			73.74%	55.33%
	Our dataset (scenario 1) 29.89 cm	10cm	91.67%	86.41%
	Our dataset (scenario 2)	31.4 cm	14.34 cm 92.06%	85.56%
	According to table 1, we can conclude: