J Koenemann

A Del Prete

Y Tassa

E Todorov

O Stasse

M Bennewitz

N Mansard

Whole-body Model-Predictive Control applied to the HRP-2 Humanoid

Controlling the robot with a permanently-updated optimal trajectory, also known as model predictive control, is the Holy Grail of whole-body motion generation. Before obtaining it, several challenges should be faced: computation cost, non-linear local minima, algorithm stability, etc. In this paper, we address the problem of applying the updated optimal control in real-time on the physical robot. In particular, we focus on the problems raised by the delays due to computation and by the differences between the real robot and the simulated model. Based on the optimal-control solver MuJoCo, we implemented a complete model-predictive controller and we applied it in real-time on the physical HRP-2 robot. It is the first time that such a whole-body model predictive controller is applied in real-time on a complex dynamic robot. Aside from the technical contributions cited above, the main contribution of this paper is to report the experimental results of this première implementation.

I. INTRODUCTION

A. Motivations

This papers deals with the problem of generating whole body motion with contacts in real time for humanoid robots. This generalized form of locomotion and manipulation would allow such robots to manipulate heavy objects, climb on uneven terrain, use the environment to stay balanced while in multiple contacts (as in Fig. 1) with one single controller. To achieve this the robot's dynamical model, forces, motor torques, self-collision, multi-contact dynamical balance have to be considered all together on a time horizon. This makes the problem very hard compared to classical instantaneous inverse kinematics. MuJoCo [START_REF] Todorov | MuJoCo: A physics engine for model-based control[END_REF], [START_REF] Tassa | Synthesis and stabilization of complex behaviors through online trajectory optimization[END_REF] realizes this by using model predictive control and a specific contact model. So far it was mostly tested in simulation [START_REF] Tassa | Control-limited differential dynamic programming[END_REF]. This paper describes a first evaluation towards a full deployment of this approach on a real humanoid robot, namely HRP-2.

B. State of the art

Whole-body motion is typically generated by inverse kinematics (IK) [START_REF] Baerlocher | Task-priority formulations for the kinematic control of highly redundant articulated structures[END_REF]. By keeping the model simple, IK is computationally efficient, but it cannot express the whole variety of constraints and tasks that we would typically expect from dynamic robots such as humanoids. For instance, the balance of a humanoid robot is an objective that is 1 Jonas Koenemann and Maren Bennewitz are with University of Freiburg jonas.koenemann@yahoo.de, maren@informatik.uni-freiburg.de 2 Andrea Del Prete, Olivier Stasse and Nicolas Mansard are with LAAS/CNRS, Toulouse adelpret@laas.fr, ostasse@laas.fr, nmansard@laas.fr 3

Yuval

Tassa is with Google UK, London yuval.tassa@gmail.com 4 Emo Todorov is with the University of Washington todorov@cs.washington.edu Fig. 1: Whole-body multi-contact experiment on HRP-2. The robot reaches for a target while it uses a table as an additional support to keep balance.

difficult to express with IK only, except by limiting the range of possible movements to only quasi-static ones. To extend the range of expressiveness of motion-generation techniques, the robotics community is switching to inverse dynamics, i.e. computing the instantaneous joint torques leading to the satisfaction of some operational constraints [START_REF] Park | Control strategies for robots in contact[END_REF], referred to as operational-space inverse dynamics (OSID) in the following. Once more, the computation cost of the method is kept low by considering only the instantaneous linearization of the system at the current time instant. This more-complex OSID model provides more expressiveness than IK: for example, it is possible to consider instantaneous balance criteria like the ZMP on flat terrain [START_REF] Vukobratović | Zero-moment pointthirty five years of its life[END_REF] or the positivity of the contact forces in case of multiple contacts [START_REF] Saab | Dynamic Whole-Body Motion Generation under Rigid Contacts and other Unilateral Constraints[END_REF]. However, the instantaneous linearization once more limits the range of expressiveness. Continuing with the same example, balance control needs to consider the future evolution of the system, because the best definition of balance might be the capability of not falling in the future [START_REF] Wieber | Viability and predictive control for safe locomotion[END_REF]. Basically, optimal control can be seen as an extension of OSID that considers the future evolution of the system [START_REF] Geoffroy | From Inverse Kinematics to Optimal Control[END_REF], rather than its instantaneous linearization.

The main difference is that the problem becomes nonlinear (and in general nonconvex) and thus much harder to solve.

In particular, it is not longer possible to guarantee that the solver outputs the global minimum of the problem (not even that any global minimum exists [START_REF] Laumond | Optimality in robot motion: optimal versus optimized motion[END_REF]): we have to accept to work with a local minimum. The properties of this local optimum highly depends on the kind of nonlinear heuristic that was used to search it. A first class of heuristic accepts a time-consuming search among unfeasible solutions to finally output a solution that might be far from the provided initial guess [START_REF] Mordatch | Discovery of complex behaviors through contact-invariant optimization[END_REF], [START_REF] Mombaur | Using optimization to create self-stable human-like running[END_REF]. In that case, the optimal-control solver rather corresponds to a path planner, providing complex solutions, but taking minutes or hours of computation. On the other hand, it is possible to always keep a feasible solution during the search; in that case, the search algorithm is realtime (we can interrupt it at any time and get a feasible trajectory) and usually much faster. The trade-off is that the optimum is often of lower quality. In that case, the optimalcontrol solver rather corresponds to a controller, which is able to track a planned trajectory, but has only limited exploration capabilities [START_REF] Tassa | Synthesis and stabilization of complex behaviors through online trajectory optimization[END_REF].

C. Approach and contributions

We focus here on a specific optimal-control solver of this second class, named differential dynamic programming (DDP) [START_REF] Jacobson | Differential dynamic programming[END_REF]. DDP is a method to exploit the sparsity of the optimal-control problem without explicit constraints. It is implemented in the software MuJoCo [START_REF] Todorov | MuJoCo: A physics engine for model-based control[END_REF]. Using multicore parallelization and a smooth dynamics approximation [START_REF] Todorov | Implicit nonlinear complementarity: A new approach to contact dynamics[END_REF], MuJoCo can compute a 500ms trajectory for a 25-DOF humanoid robot in 50ms on a standard desktop computer. Based on this efficiency, it is possible to update in real-time the robot trajectory and control it from the output of the nonlinear optimal-control solver. This is referred to as model predictive control (MPC) [START_REF] Alamir | Stabilization of nonlinear systems using receding-horizon control schemes: a parametrized approach for fast systems[END_REF]. As described in [START_REF]Fast NMPC: A reality-steered paradigm: Key properties of fast NMPC algorithms[END_REF], MPC proposes some additional challenges with respect to standard optimal control. In particular, MPC has to be able to handle the noise due to the robot's sensors and the delay due to the computation time.

The contribution of this paper is to report the first application of whole-body nonlinear model predictive control on a humanoid robot (i.e. HRP-2). It is the first time that MPC is applied to the whole body of a humanoid robot, controlling the balance, handling different constraints (e.g. collisions, joint limits) and performing the specified reaching task at the same time.

In the next Section II, the MPC scheme and the DDP solver are introduced. In particular we describe choices and technical contributions that we developed to apply DDP as an MPC on a physical robot. In Section III, the connection between torque-based MPC model and physical positiondriven actuators is described. The experimental setup is detailed in Section IV and the results are described and commented in Section V. Section VI draws the conclusions.

II. MODEL PREDICTIVE CONTROL

A. Challenges in model predictive control

The optimal control problem (OCP) consists in finding the control and state trajectories optimizing a cost functional, typically composed of an integral term and a terminal term. OCP outputs an optimal state trajectory X * that can be followed in open-loop from the starting point to the goal.

Of course, open-loop tracking raises significant issues as soon as the predicted model does not perfectly match the physical system. The general idea of MPC is to update the OCP at every iteration of the control loop and only apply on the physical system the controls corresponding to the very first part of the trajectory [START_REF] Diehl | A real-time iteration scheme for nonlinear optimization in optimal feedback control[END_REF]. MPC is a key methodology to properly control complex robots such as humanoids. It would for example unlock behaviors such as multi-contact locomotion.

The OCP is a planning problem, with its own difficulties of existence and approximation of the solution. Besides these problems, MPC raises additional specific challenges. First, in general we can not neglect computation times, and more generally neglect the delays between the sensor measurements and the computation of the associated control [START_REF] Zavala | The advanced-step nmpc controller: Optimality, stability and robustness[END_REF], [START_REF] Alamir | Monitoring control updating period in fast gradient based nmpc[END_REF]. We will see later that the empirical delay obtained on HRP-2 corresponds to 10 measurement-and-control cycles of the real robot. The use of the sensor measurements should be carefully considered with respect to this asynchronous delay. Several measurements can be merged to predict the robot state before each new computation. The state prediction must take into account the current set of optimal control candidates [START_REF]Fast nmpc: Some good news and some facts to keep in mind[END_REF]. Additionally, each measurement can also be immediately used at the next control cycle, to servo the robot on the lastly-computed optimal trajectory. The DDP is particularly interesting for this last issue, because it outputs the optimal Riccati gains in addition to the optimal trajectory [START_REF] Li | Iterative linear-quadratic regulator design for nonlinear biological movement systems[END_REF], providing an locally-optimal policy to feedback on new sensor measurements while waiting for the DDP solver to output a new trajectory candidate.

The control literature already described these three challenges and proposed some solutions, but the corresponding questions are still open. The main contribution of this paper is to experimentally exhibit these problems for the first time on a humanoid robot. In this section, we first recall the main facts about DDP that are necessary to understand the rest of the paper. We then describe with more details the three relevant problems that we noticed when applying MPC to the real robot, and quickly describe the preliminary solutions that we have applied in the experiments.

B. Differential dynamic programming

We consider in the following a direct resolution of OCP, i.e. iteratively approximating a local optimum to the OCP by a sampled (discrete-time) trajectory of state and control. This direct method boils down to the resolution of a large nonlinear optimization problem, using typically a Newton or quasi-Newton method. One of the important feature to take into account in the implementation of such a direct OCP solver is the intrinsic sparsity arising from the temporal structure of the problem. DDP [START_REF] Jacobson | Differential dynamic programming[END_REF] is an efficient way to handle this sparsity while keeping the implementation simple, in the particular case where the OCP is unconstrained. In that case, it is nearly equivalent to a Newton descent [START_REF] Murray | Differential dynamic programming and Newton's method for discrete optimal control problems[END_REF]. We considered a quasi-DDP descent named iterative LQR (iLQR), arising when neglecting the secondorder derivatives of both the dynamics function and the cost residual [START_REF] Li | Iterative linear-quadratic regulator design for nonlinear biological movement systems[END_REF]. In essence, this is very similar to the Gauss-Newton optimization descent. iLQR is particularly appropriate for control, since the second-order approximation increases the computational efficiency (with the drawback, not very important for control, that the exploration capabilities are lowered by the lack of second-order terms). iLQR has largely been described in previous papers [START_REF] Li | Iterative linear-quadratic regulator design for nonlinear biological movement systems[END_REF], [START_REF] Tassa | Synthesis and stabilization of complex behaviors through online trajectory optimization[END_REF], [START_REF] Tassa | Control-limited differential dynamic programming[END_REF]. We recall the major technical facts in the Appendix. Mainly, the reader should be aware that DDP: i) is a Newton-like optimization scheme using only first-order derivatives, ii) takes advantage of the sparsity of the derivatives to keep a reasonable computation cost, and iii) is able to consider box constraints on the control variables. The solver outputs the optimal state trajectory X * and the corresponding control U * , but also the gains K of the optimal LQR to regulate X. We call Π = (U, K) the output policy.

C. Cost overview and real-time property

The cost of one iteration of the DDP solver is in O(m 3 N ), with m being the dimension of state and control, and N being the number of preview steps. In the experiments, this typically represents 10ms of computation (see Section VI for details). The solver additionally needs to compute the derivatives of both the dynamics and the cost function. For the humanoid robot in contact with the environment, the dynamics are computed by a contact solver and have no closed-form; the derivatives can only be approximated by finite differences. Thanks to the parallel computation on several cores, this typically takes 40ms in the experiments. These two parts are real-time: the computation time is bounded and typically has low variance. On the other hand, the solver also performs a line search (see Appendix), for which we do not know in advance the number of iterations. Each iteration is in O(N ) and typically takes around 1ms in the experiments. The intermediate solutions of the line search are not valid because they cannot be applied on the robot: this final part of the algorithm is therefore not real-time.

In summary, despite the parallelization, the computation time of one iLQR iteration (around 50ms in the experiments) is not negligible with respect to the robot control rate (typically around 5ms) and it is not real-time. However, it is fast enough to be recomputed online (driving the motion of the robot) in a second thread running independently of the real-time control thread.

D. Real-time and fast-time loops

The robot control loop is typically real-time (the motors inputs must be updated at a regular and fast rate). On the opposite, the optimal-control computation is typically slower and cannot be guaranteed to converge in real-time, i.e. within a fixed time. The optimal-control solver must then be embedded in a second thread, which is independent of the real-time control thread. We say that this second thread is fast-time, i.e. it cycles as fast as possible, with some nonconstant rate variance. It is not necessary that this second thread runs on the same computer as the control thread (in our case, we used a laptop-style CPU onboard the robot to host the control thread, and a desktop-style CPU offboard for the optimal-control solver).

The sensor measurements are received by the real-time control thread, and transferred to the optimal-control thread. Both the trajectory optimization and the network communication introduce a delay τ , named the policy lag. The lag MuJoCo HRP-2

t i t i+1 t i+2 x(t i ) x(t i+ 1 ) x(t i+2) ⌧ i ⌧ i+1
Ui( ti

: ti + T ) Ui 1 ( ti 1 : ti 1 + T ) Ui+1 ( ti+1 : ti+1 + T ) U i 1 ( ti : t i+1 ) U i (t i+1 : ti+1 ) ti ti 1 ti = t i + ⌧ ti+1
Fig. 2: Example of communication between the real-time control thread ("HRP-2") and the DDP thread ("MuJoCo"). At time t i , the updated policy Π i-1 is received before expected (τ i-1 < τ ): the application of the new policy is delayed. At time t i+1 it is the opposite (τ i > ): after discarding the first outdated samples, the updated policy is immediately applied. Again, at time t i+2 , the lag is overestimated, so the application of the updated policy is delayed.

has a varying length, depending on hardware and algorithm behavior.

E. State prediction

Each cycle of the optimal-control thread is denoted by the index i. It corresponds to a measurement of the robot xi captured at time t i , and to an output policy Π i . The policy is received by the control thread after the lag τ i . Since t i + τ i corresponds to the time of the next-cycle measurement, we denote it by t i+1 . The policy Π i can be applied at best from t i+1 . It is therefore not necessary to try to optimize the controls Π i during the policy lag (i.e. from t i to t i+1 ) since the robot will not be able to apply them. By using a proper estimation τ of τ i , we therefore predict the state at t i + τ by integrating the previous policy Π i-1 starting from the measurement xi .

When τ is bigger than the actual τ i , the policy switch is delayed (when the policy arrives at time t i+1 , the corresponding control has not been optimized and we have to wait until t i + τ to have an updated control). When τ is smaller than the actual τ i , the policy switch immediately applies, but the first controls are not used since they are outdated. A schematic view of the policy lag and state prediction is given in Fig. 2.

F. Tracking the state

The DDP outputs redundant information: the control policy (3) (composed of the feedforward term k and the feedback matrix K) and the corresponding state trajectory x * . The control to send to the motor should then be:

u = -k -K(x -x * )
, where xx * is the deviation from the planned trajectory. The feedback gain K will then "distribute" the correction, depending on the effect of the deviation on the cost. Alternatively, it is possible to neglect the feedback and only apply the feedforward u = -k. Or only the reference state x * can be used as input to a separate controller. This last approximation has the advantage to lower the transmission bandwidth (K and k are not transferred). It is also interesting when the robot has a very good tracking system, or when the low-level control loop is difficult to model (and therefore to integrate in the dynamics of the MPC). This last point is discussed in Section IV.

We experimentally tested several ways of tracking the optimum computed by the MPC and we report the conclusions in the experimental Section VI.

G. MPC overview

Finally, the two main loops on the real-time control thread and on the optimal-control threads are summarized respectively in Alg. 1 and 2. In addition to these two loops, a "measurement" server runs on the optimal-control thread to receive the measurements sent from the control thread. The measurements are buffered when they arrive but most of them are simply discarded: only most recent measurement before starting a new optimal-control cycle is used. And a "policy" server runs on the real-time control thread to receive the optimum computed by the DDP. At the end of the transition interval, the second-last policy is discarded. Typically the forward dynamic model of a mechanical system takes forces as inputs. This implies that we should be able to control the joint torques on our humanoid robot. However, HRP-2 has high-ratio gear boxes, which introduce large friction (especially static friction), making torque control challenging [START_REF] Traversaro | Inertial Parameter Identification Including Friction and Motor Dynamics[END_REF]. Implementing torque control on this kind of robots requires a torque feedback, which is currently not available on HRP-2. That is why HRP-2 is equipped with a high-gain position control, which allows for accurate tracking of the desired joint trajectories, but makes the robot stiff.

1 for i = 0 till end,
To account for this mismatch, we modeled the position controller in the forward dynamics of the system. On each motor a fast control loop ensures the tracking of the desired motor current i d , which is computed as:

i d = K p (u -q) -K d q,
where q is the joint position and K p , K d are the proportional and derivative gains, respectively. The control input u is then the desired joint position commanded to the low-level controller. The motor current i is directly proportional to the motor torque µ m , which -neglecting the gear-box elasticity -differs from the joint torque µ (mainly) because of the friction torque µ f :

µ = K m i µm -K v q -µ s µ f
, where we split the friction torque into two parts: viscous friction K v q and static friction µ s . Under the assumption of perfect current tracking we have then:

µ = K m K p K1 (u -q) -(K m K d + K v ) K2 q -µ s
To add the actuator dynamics inside the robot model we would need to identify K 1 , K 2 and µ s . Unfortunately modeling static friction is known to be a challenging problem [START_REF] Olsson | Friction models and friction compensation[END_REF] because it depends on position, temperature and external forces. However, due to the high gains, the position control compensates for most of the static friction, justifying our choice to neglect it. The identification is performed by exciting the system with sinusoidal trajectories on each joint, while collecting the desired positions u and real positions q meas . Then a nonlinear least-squares problem is solved to identify the gains K 1 , K 2 by minimizing the squared difference between the measured trajectory x meas = q meas q meas and the one obtained by simulating the system:

min K1,K2 ||x -x meas || 2 s. t. x i+1 = f (x i , u i , K 1 , K 2 )
For the simulation of the system we used the MuJoCo simulator, whereas for the optimization we used the Matlab function lsqnonlin.

The choice of modeling the position controller and the actuator inside the dynamics of our system has a number of pros and cons. On the bright side, MuJoCo is inverting the actuator-control model, simplifying the work of the onboard controller. Moreover, by bounding the control inputs we can ensure the respect of the joint limits. Furthermore, we could use a model-based reinforcement learning approach to identify on the fly the parameters of our model, while generating the optimal control policy. On the other side, in this case the actuator-control dynamics is rather simple and we could also invert it outside MuJoCo. Also, bounding the control inputs inside the joint limits could prevent the robot from applying the desired forces on the environment. Finally, the full range of consequences of this choice are still unclear to us and we believe they need further investigation.

IV. EXPERIMENTAL SETUP

A. Hardware

We carried out all the tests on the humanoid robot HRP-2.

In the first tests we used only the upper body of the robot (10 DOFs, 4 in each arm and 2 in the torso). In the last tests we used instead its whole body (27 DOFs), thus dealing with an underactuated system. Onboard the robot we run a realtime linux on a 2.93 GHz CPU where the real-time position controller loops every 5ms, using only encoder feedback. The local stability of the robot is enforced by the proprietary stabilizer, which exploits the IMU and force/torque-sensor feedback. The machine used for the MPC computation is offboard the robot and it runs a non-real-time Windows on a 12-core 4 GHz CPU. The two machines are connected through a standard Wi-Fi (round-trip time below 1ms).

B. Software

In all our tests we used MuJoCo [START_REF] Todorov | MuJoCo: A physics engine for model-based control[END_REF] for the online trajectory optimization. Inside MuJoCo a server thread waits for receiving the state from the robot, while a client thread sends the computed trajectories to the robot. Onboard the robot, a non-real-time server thread waits for the updated control sequence, takes care of the interpolation (MuJoCo's trajectories have a 20ms time step, while the control period is 5ms) and it updates the buffer containing the future sequence. At the same time, the real-time loop takes care of the PD joint-position control, which applies the last control received from MuJoCo.

C. Overview of the Experiments

In the first test the robot had to use its upper body to reach a moving target. The target (i.e. a yellow ball) existed in simulation only and it was moved by the user through a 3D mouse. Using only the upper body we could neglect the balancing problem, and focus on validating the position control and the state prediction. Then in the second test the robot used its whole body to reach the same moving target while balancing (see Fig. 6). In the third test we ask the robot to make contact on a table with its right hand, and then to reach for a ball in front of it with the left hand (see Fig. 1).

V. RESULTS

We tested the presented control architecture in different scenarios, both in simulation and on the real robot. The simulation results led us to modify the MPC scheme to overcome two issues: discontinuities in the control sequence and divergence of the DDP algorithm. Only the modified MPC scheme was then tested on the real robot.

A. Position-Driven Model -Upper-Body Reaching

This test validates the high-gain position control that we inserted inside the dynamical model of the system. Fig. 3 shows the values (scaled to properly fit the plot) of the cost terms along time, which are: i) distance between hand and target, ii) z axis of torso (to keep it vertical), iii) joint velocities, iv) joint torques. Additionally, the whole-body dynamics and the contact constraints are enforced by the optimal-control solver, as well as the collision and joint-limit contraints.

The experiment on the real robot (synchronized with the movement planned by the solver MuJoCo) is presented in the accompanying video 1) State Update and State Prediction: The first problem that we faced have been the discontinuities in the control sequence introduced by the state update. We observed a discontinuity every time the controller switched between two subsequent control sequences received from MuJoCo. Fig. 4 shows that these discontinuities are strongly attenuated when not updating the state. This led us not to use any state update in the experiments on the robot. To attenuate the discontinuities even more we performed a state prediction (see Section II-E). In other words, each optimization starts from the state predicted by the integration of the dynamics (considering the policy lag) rather than from the current state. Fig. 5 shows that the state prediction that we implemented actually helps reducing the discontinuities of the control sequences. Of course this approach would not work in case of strong external disturbances, but it allowed us to eliminate the above-mentioned discontinuities and focus on other crucial aspects.

B. Torque-Driven Model -Whole-Body Reaching

The task was expressed through a cost function that trades off the following (conflicting) objectives: i) reach the target with the right hand, ii) keep the feet at their initial position, iii) keep the capture point at the center of support polygon, iv) keep the orientation of torso/waist upright, v) penalize joint velocities, iv) penalize joint torques. Like previously, the dynamics, contact and collision constraints are also enforced by the solver. Fig. 7 shows the values of some cost terms, while Fig. 6 shows some screenshots of both the MuJoCo simulation and the real robot. The whole sequence can be watched in the accompanying video. 1) Actuator Model: In this test we could not use the position-driven model in MuJoCo because it led to divergence of the DDP. We believe that this did not happen in the first test because of the lack of contacts with the ground. To overcome this issue we used a standard torquedriven dynamic model; then we sent the state trajectories as references for the joint-position control.

2) Self-Collision Avoidance: Fig. 8 proves the effectiveness of the self-collision avoidance implemented in MuJoCo. Even if we moved the target inside the waist of the robot we see no penetration between the capsules used to approximate the robot's links.

C. Torque-Driven Model -Whole-Body Multi-Contact

Finally, we performed a reaching motion while controlling the robot balance with the other hand. This experiment confirms the capabilities of the optimal-control approach, that is able to unify many different operation modes of the robot in a single controller. For this test the cost function of the previous test has been marginally modified to bring the robot in contact. We removed the term penalizing inclinations of the torso/waist and we decreased the weight of the capturepoint term. Moreover, we added a term to reach a target with the left hand. Contrary to the previous tests, here we used a sequence of four different cost functions to guide the motion of the robot over time. The difference between the cost functions was only the target positions of the right/left hands (see position of the colored balls in Fig. 1). To get the right hand in contact with the table we moved the associated target below the table. The whole sequence can be seen in the accompanying video.

VI. CONCLUSIONS

This paper described the first application of real-time whole-body MPC on a full humanoid robot. We carried out all the experiments on the HRP-2 robot using MuJoCo for the online trajectory optimization. Besides the experimental results, the contributions of this work are the modifications of the MPC scheme that we implemented to address two practical issues: the discontinuities in the control sequences and the actuation differences between model and real robot. These issues appear only when working with a real complex mechanical system, which is why they have been overlooked by previous simulation studies.

In the near future we will investigate how to reintroduce the state update inside the MPC scheme while maintaining smooth control sequences. Also we are in the process of implementing joint-torque control on HRP-2, which would allow us to use a torque-driven model for the trajectory optimization. Our goal is to gradually increase the complexity of the generated motion: our next step is to have the robot making dynamic motion while making and breaking contacts with the environment, e.g. walking with additional support of the hands.

APPENDIX DIFFERENTIAL DYNAMIC PROGRAMMING A. Problem definition

Consider a discrete dynamical system with state x and control u:

x i+1 = f (x i , u i ) (1)
The optimal control problem consists of minimizing a userdefined cost J(x 0 , U ) over a certain time horizon N :

U * = argmin U N -1 i=0 l(x i , u i ) + l f (x N ) J(x0,U )
Since DDP is an implicit method, it represents the trajectory through the control sequence U = {u 0 , . . . , u N -1 }, while the states are computed integrating [START_REF] Todorov | MuJoCo: A physics engine for model-based control[END_REF]. Let us define the cost-to-go as:

J j (x j , U j ) = N -1 i=j l(x i , u i ) + l f (x N ),
where U j = {u j , . . . , u N -1 }. The Value at time i is the optimal cost-to-go starting at x i :

V i (x i ) = min Ui J i (x i , U i ) = min ui [l(x i , u i ) + V i+1 (f (x i , u i ))]

B. Dynamic Programming

The Dynamic-Programming principle simplifies this minimization over the entire control sequence U to a cascade of minimizations over the single controls u i . Starting by setting V N (x N ) = l f (x N ), we can optimize backwards in time until we find the entire optimal control sequence U * .

C. Differential resolution

Solving this cascade is not possible in the general case, due to the well-documented Bellman's curse of dimensionality. To solve each nonlinear optimization, the cost variation is therefore approximated with a quadratic function: Q(δx, δu) =l(x, u) -l(x + δx, u + δu)+ V (f (x, u)) -V (f (x + δx, u + δu))

≈   1 δx δu     0 Q x Q u Q x Q xx Q xu Q u Q ux Q uu     1 δx δu   , (2)
where the coefficients (i.e. first and second derivatives of Q with respect to x and u) are simply obtained as functions of the derivatives of the cost l(x, u), the Value V (x) and the dynamics f (x, u) (see [START_REF] Tassa | Synthesis and stabilization of complex behaviors through online trajectory optimization[END_REF] for the exact expressions). While DDP uses the complete expressions of these coefficients, iLQG neglects the second derivatives of the dynamics in order to make the computation faster. In both cases, the solution of the quadratic approximation is then:

δu * = argmin δu Q(δx, δu) = -Q -1 uu Q u k -Q -1 uu Q ux K δx (3)
Substituting δu * in (2) we get a quadratic model of V , which we need to propagate backwards in time.

D. Line search

After computing δu * for the whole horizon we update the control policy Û and the corresponding state trajectory X starting from x0 = x 0 :

ûi = u i -αk i -K i (x i -x i ) xi+1 = f (x i , ûi ),
where α is a backtracking search parameter, initially set to 1 and then iteratively reduced. The complete description of the nonlinear heuristic to adjust α is based on the Levenberg-Marquardt algorithm and it is described in [START_REF] Tassa | Synthesis and stabilization of complex behaviors through online trajectory optimization[END_REF].

Fig. 3 :

 3 Fig. 3: Cost terms during upper-body reaching.

  Control sequence without state update.

Fig. 4 :

 4 Fig. 4: Control sequences for the 4 joints of the left arm during a reaching motion. The discontinuities are much larger when updating the state.

  With state prediction, expected policy lag τ = 0.12 s. Without state prediction, expected policy lag τ = 0 s.

Fig. 5 :

 5 Fig. 5: Numerical derivative of the control sequences for the 12 leg joints during a reaching motion. Derivatives are computed by finite differences. Blue circles represent the difference between the policy lag and the expected policy lag. When this is greater than zero a discontinuity is expected.

Fig. 6 :

 6 Fig.6: Whole-body experiment on the real robot. The robot tracks a ball with its right hand while it keeps balance. The whole sequence can be seen in the accompanying video.

Fig. 7 :

 7 Fig. 7: Cost terms during whole-body reaching.

Fig. 8 :

 8 Fig. 8: Distances between the closest capsules. Negative distances would imply penetration.


None Bertrand Bardet

None; N Anabelle Sequeira-Le Grand

Niyazi Acar

Catherine P Garcher

Alain M Bron

Lionel Brétillon

L Bretillon

None; C P Acar

None A M Garcher

None; L Bron

None. Bretillon

Retinal ganglion cell loss evaluation by flow cytometry after intravitreal injection of endothelin-1 in the rat

Keywords:

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


Pierre Buyssens

Isabelle Gardin

Su Ruan

Abderrahim Elmoataz

Eikonal-based Region Growing for Ecient Clustering

Keywords: Superpixels, segmentation, clustering, Eikonal equation

We describe an Eikonal-based algorithm for computing dense oversegmentation of an image, often called superpixels. This oversegmentation respect local image boundaries while limiting undersegmentation. The proposed algorithm relies on a region growing scheme, where the potential map used is not xed and evolves during the diusion. Renement steps are also proposed to enhance at low cost the rst oversegmentation. Quantitative comparisons on the Berkeley dataset show good performance on traditional metrics over current state-of-the art superpixel methods.

Introduction

With the increasing amount of available data, and the need for fast and accurate processing, the simplication of data becomes a crucial point for many applications. A convenient way to address this task is to consider that datas can be modeled with a graph G = (V, E, w), where V is the set of vertices, E a set of edges, and w > 0 is the weight function that models the interaction between vertices. Exhibiting clusters of this graph leads to a simplication of the data and decreases the size of the problem.

Many techniques of graph clustering have been proposed such as cut-based, spectral or random walk methods (see [START_REF] Satu | Graph clustering[END_REF] for a comprehensive review of these techniques).

Recent works [START_REF] Desquesnes | Eikonal equation adaptation on weighted graphs: Fast geometric diusion process for local and non-local image and data processing[END_REF] adapt the eikonal equation to graphs in order to perform over-clustering from an initial set of annotated vertices V 0 . Let f : V → R be a real-valued function that assigns a real value f (u) to each vertex u ∈ V .

The reformulation of the Eikonal equation in the graph domain leads to the equation :

   (∇ f w )(u) p = P (u) ∀u ∈ V f (u) = φ(u) ∀u ∈ V 0 (1)
where (∇ - w f )(u) is the weighted morphological gradient at a vertex u (see [START_REF] Lézoray | Image Processing and Analysis with Graphs: Theory and Practice[END_REF] for details), P is a positive function, and φ is an initialization function.

In this paper, we focus on grid graphs for image processing with the aim of grouping perceptually and adjacent pixels into meaningful regions, the so-called superpixels. Superpixels have become an important step in many computer vision applications such as segmentation [START_REF] Gould | Multi-class segmentation with relative location prior[END_REF][START_REF] Yang | Layered object detection for multi-class segmentation[END_REF], object localization [START_REF] Fulkerson | Class segmentation and object localization with superpixel neighborhoods[END_REF], depth estimation [START_REF] Zitnick | Stereo for image-based rendering using image over-segmentation[END_REF], and scene labeling [START_REF] Farabet | Learning hierarchical features for scene labeling[END_REF].

Some properties of an algorithm that generate superpixels are often desirable : (1) Superpixels should adhere well to object boundaries while limiting undersegmentation errors, [START_REF] Arbelaez | Contour detection and hierarchical image segmentation[END_REF] as superpixel methods are used as preprocessing, the algorithm should have a low complexity, (3) it has to be simple to use (i.e. few parameters). In addition, some other properties may be desired : the control of the amount of superpixels, or the compactness of them.

Several superpixel algorithms exist in the literature, they can be roughly divided into two approaches : The rst consists in gradually growing superpixels from an initial set of centers. This approach includes Watershed [START_REF] Vincent | Watersheds in digital spaces: an ecient algorithm based on immersion simulations[END_REF],

Turbopixels [START_REF] Levinshtein | Turbopixels: Fast superpixels using geometric ows[END_REF], SLIC [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF], Consistent Segmentation [START_REF] Cw Zitnick | Consistent segmentation for optical ow estimation[END_REF] and Quick Shift [START_REF] Vedaldi | Quick shift and kernel methods for mode seeking[END_REF]. The second approach relies on a graph formulation of the problem and aims at nding an optimal cut according to an objective function that takes similarities of neighboring pixels into account. This approach includes Entropybased energy function method [START_REF] Liu | Entropy rate superpixel segmentation[END_REF], optimal cuts [START_REF] Moore | Superpixel lattices[END_REF][START_REF] Moore | lattice cut'-constructing superpixels using layer constraints[END_REF], graph-cut [START_REF] Veksler | Superpixels and supervoxels in an energy optimization framework[END_REF], and agglomerative clustering of the nodes of the graph [START_REF] Felzenszwalb | Ecient graph-based image segmentation[END_REF].

In this paper, we propose a new algorithm for superpixel generation :

Eikonal-based Region Growing Clustering 1 (ERGC) that starts from an initial set of seeds and dilates them, and then renes the result oversegmentation by adding/moving cuts. It formulates the superpixel segmentation task as a solution of an Eikonal equation. Equation 1 becomes :

   ∇U (x) = F (x) ∀x ∈ I U (x) = 0 ∀x ∈ Γ ( 2
)
where I is the image domain, F a positive function, Γ the set of initial seeds, and U (x) the traveling time or geodesic distance of x from source Γ. Focusing on grid graphs, it can be solved eciently with the Fast-Marching method.

The major change proposed in this paper concerns the function F , which is not xed and evolves during the front evolution. It is detailed at Section 2.2.

1 Source code and executable of ERGC can be found at https://sites.google.com/site/pierrebuyssens/ergc ERGC is simple to use (by default, the only parameter is the desired number of superpixels), as fast as other superpixels methods, and outperforms them on two of the three traditional metrics.

The rest of the paper is organized as follows : Section 2 details the proposed potential function F , and the ERGC algorithm. Section 3 gives qualitative and quantitative comparisons of performances between ERGC and state-of-the-art methods. Some aspects and extensions of the proposed method are than discussed in Section 4, while Section 5 concludes the paper.

Superpixels method 2.1. Notations

In the following we adopt several notations to simplify the reading of the paper. A particular pixel of image I is noted p and consists of a coordinate couple (x p , y p ). A region R i consists of a seed pixel s i and a size N i in pixels.

The color of a pixel p is noted C p , and the mean color of a region R i is noted

C i .
Note that the color images are considered in the CIELAB colorspace, so the color vector of a pixel (or a region) C reduces to [l, a, b] T .

Proposed potential function

Since a superpixels method aims at grouping perceptually and adjacent pixels into meaningful regions, we propose a potential function F that conveys this desirable property. The right term of equation 2 is computed according to the mean color of the adjoining region :

F c (p, R i ) = C p -C i 2 2 (3)
This potential function measures the perceptual color distance between the pixel p and the region R i . For color images in the CIELAB colorspace, F reduces to :

F c (p, R i ) = (l p -l i ) 2 + (a p -a i ) 2 + (b p -b i ) 2 (4)
where [l i , a i , b i ] T is the mean color vector of region R i .

In comparison to traditional gradientbased approaches [START_REF] Desquesnes | Eikonal equation adaptation on weighted graphs: Fast geometric diusion process for local and non-local image and data processing[END_REF] where F (p) = ∇I , the proposed formulation favors the grouping of similar pixels, even for pixels that are far from the initial seeds (Fig. 2).

As a numerical solver of the Eikonal equation 2, we adopt the Fast Marching method introduced by Sethian in [START_REF] Sethian | Level set methods and fast marching methods: evolving interfaces in computational geometry[END_REF]. It uses a priority queue to order the pixels as being the current estimate of the geodesic distance to the closest seed (see [START_REF] Peyré | Geodesic methods in computer vision and graphics[END_REF] for a detailed description of the fast marching algorithm).

Within the fast marching algorithm, each time a pixel p is inserted to a region R i , the attributes of this region are easily updated :

   C i ← C i ×N i +Cp N i +1 N i ← N i + 1 (5)
This formulation clearly makes the potential of a pixel dependant on the features of an adjoining region, which is updated during the process. It promotes the diusion to pixels whose color is close to the color of the region, hence creating homogeneous superpixels. A region can also absorb smoothly localized noisy pixels since such a pixel contributes weakly to the mean color of the region.

U (p) =    1 2 (U (q) + U (r) + √ ∆) if ∆ ≥ 0 min(U (q), U (r)) + F otherwise ( 6
)
where F is computed with equation 3, q and r are the neighbors of p such that

U (q) = min(U (x -1, y), U (x + 1, y)), U (r) = min(U (x, y -1), U (x, y + 1)),
and ∆ corresponds to solving the equation (u

-U (r)) 2 + (u -U (q)) 2 = F 2
(see [START_REF] Peyré | Geodesic methods in computer vision and graphics[END_REF] for details). if p is a seed then 3:

Σ p ← ALIVE, U (p) = 0, add p to L 4: else 5: Σ p ← FAR AWAY, U (p) = ∞ 6:
end if 7: end for 8: while L = ∅ do 9:

p ← arg min q (U (q) | q ∈ L) 10:

Σ p ← COMPUTED 11:
Update adjoining superpixel R i (equation 5)

12:

for q ∈ N eigh(p) do 13:

compute F (q, R i ) (equation 3)

14:

compute local solution t(q) (equation 6)

15:

if Σ q = ALIVE and t(q) < U (q) then 16:

U (q) = t(q) 17: end if 18: if Σ q = FAR AWAY then 19:
Σ q ← ALIVE 20:

U (q) = t(q) 21: Add q to L 22: end if 23:
end for 24: end while decreases as the potential is high. Nevertheless, after a given time, the front ends up passing through the contour, and the potential becomes low again, so the front can evolve with a high speed. This behavior is not a surprise and can be seen on the red-blue image of Figure 2. The proposed potential function adds to the diusion a sort of memory of the initial color of the seeds.

When a front arrives on a region with a dierent color, the potential becomes very high. Even if the front passes through the interface separating the two regions, the potential remains very high, hence prohibiting the diusion too much.

Algorithm initialization

The initialization consists of sampling K seed pixels on a regular grid with an interval S equal to S = N/K with N the number of pixels in the image (initialization similar to [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF] and [START_REF] Levinshtein | Turbopixels: Fast superpixels using geometric ows[END_REF]).

The mean and variance color of the seed pixel and its 4-connexity neighboring pixels is then computed. The same computation is performed for pixels that lie in a 3 × 3 neighborhood, and the seed is moved to the pixel that lowers the variance color. Such a perturbation of the initial seeds avoids potential outlier pixels as seeds, favorably initializes the diusion, and gives a more robust initial mean color for each superpixel.

Complexity and segmentation speed

After the initialization, a front propagation is performed with an online update of the superpixels. Complexity of the diusion relies essentially on the complexity of the fast marching algorithm, which is roughly in O(n log(n))

with an appropriate heap for sorting the pixels according to their geodesic distance. Despite this theoretical complexity, the proposed algorithm is very fast in practice, and is nearly linear in time. Note that using dierent data structures (and additional storage), dierent O(n) implementations have been proposed in [START_REF] Yatziv | O(n) implementation of the fast marching algorithm[END_REF][START_REF] Kim | An O(n) level set method for eikonal equations[END_REF].

Renement by adding new superpixels

In this section, we propose a renement of the clustering that iteratively adds new seeds after a full pass of the algorithm. Since initial seeds are placed on a grid, some objects of the image may not contain an appropriate seed, and these objects may not be nally well segmented (Fig. 3(c)). In such a case the resulting geodesic distance map U exhibits high values in this area (Fig. 3(d)). The renement consists of adding a new seed to the location of the maximum geodesic distance, and to recompute the solution of equation 2. Since it is unnecessary to apply the algorithm to the whole image, only a small part of the image around the new seed is considered (Fig. 3(b)). The renement is summarized as follows :

1. Perform a full pass of ERGC with the seeds placed on a grid,

Renement by moving superpixels

In this section, we propose a simple procedure to increase the quality of the oversegmentation without increasing the total number of superpixels. It consists in adding a new superpixel while removing a weak one, such that the global number of superpixels remains constant.

The selection of the superpixel to remove consists of 3 steps :

1. Given an oversegmentation of the image, compute the underlying Region Adjacency Graph (RAG), 2. For each vertex v of the graph, compute its normalized volume vol:

vol(v) = u∼v w(u, v) Card(N v ) (7)
where u ∼ v denotes two adjacent vertices, w(u, v) the weight of the edge connecting u and v, and N v the set of neighbors of v.

3. Select the superpixel corresponding to the vertex with the minimal normalized volume as the superpixel to remove.

This procedure selects the superpixel that is the closest to its neighbors.

The weight function is the L 2 norm and reects the dierence between two adjacent superpixels.

The initial oversegmentation is then rened by adding a new superpixel with the procedure 2.5 detailed above, while the superpixel to remove is simply discarded.

The whole procedure is iterated until a criterion is reached. In the following, it is stopped when I U no longer decreases. We also add a limit of 10 iterations, which in practice, is rarely reached. (from green boxes), and added at other locations (black boxes).

Figure 4 shows a renement of an initial oversegmentation by moving 3 superpixels. Both oversegmentations contain 150 superpixels. The complexity of this renement is low since it is similar to the renement by adding new seeds.

Summary of the algorithm

Although a spatial constraint can easily be added to F (see gure 3), by default the only parameter of ERGC is the number of desired superpixels.

The whole algorithm consists of 3 steps :

1. Proceed to the initial diusion with 90% of the desired seeds placed on a grid, 2. Rene the oversegmentation by adding the remaining 10% of seeds (Section 2.5),

3. Iterate the moving seeds procedure (Section 2.6) until the stopping criterion is reached.

For example, the algorithm produces 200 superpixels by rst placing 180 seeds on the initial grid, then by adding 20 more seeds with the procedure described in Section 2.5.

These ratios ensure that the whole algorithm runs in a reasonable amount of time (about half a second on a Berkeley image).

Comparison with StateoftheArt

We compare ERGC to stateoftheart methods SLIC 2 [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF], Entropy Rate Superpixels 3 [START_REF] Liu | Entropy rate superpixel segmentation[END_REF] (ERS), SEEDS 4 [18], TurboPixels 5 [10] (TP), and Gradient based diusion [START_REF] Desquesnes | Eikonal equation adaptation on weighted graphs: Fast geometric diusion process for local and non-local image and data processing[END_REF] (GrB). Examples of superpixel segmentations produced by each method appear in Figure 7.

The popular SLIC method proposed in [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF] adapts kmeans in the image plane to iteratively exhibit superpixels. The addition of a spatial constraint term produces regular regions that adhere quite well to image. The ERS algorithm proposed in [START_REF] Liu | Entropy rate superpixel segmentation[END_REF] oversegments an image via an obective function 2).

The Berkeley dataset [START_REF] Arbelaez | Contour detection and hierarchical image segmentation[END_REF] is used as a benchmark. It consists of 500 images of size 481×321 (or 321×481) and several ground truth manual segmentations for each image.

Figure 6 shows quantitative results on standard metrics, including Boundary Recall, Undersegmentation Error and Achievable Segmentation Accuracy.

We also add a compactness metric that reects the compacity of the superpixels. As a baseline for Boundary Recall, Undersegmentation Error and Achievable Segmentation Accuracy metrics, we also show the performances of a grid of square superpixels (GRID).

All results are computed from scratch using these evaluation metrics and the same hardware. Default parameters are used for the stateoftheart methods.

Boundary recall

Boundary recall (BR) measures the fraction of ground truth edges that is also present in superpixel segmentation within a distance threshold t. In our experiments, t is xed to 2 as in [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF][START_REF] Liu | Entropy rate superpixel segmentation[END_REF][START_REF] Van Den Bergh | Seeds: superpixels extracted via energy-driven sampling[END_REF]. Fig. 6(a) shows boundary recall measures for each method according to the number of superpixels. As the number of superpixels increases, the boundary recall is naturally higher.

SEEDS outperforms all other algorithms on this metric. Nevertheless, this result has to be appreciated in the light of the superpixels compactness. As shown at Fig. 6(d) (second and third rows), ERS and SEEDS superpixels have the lowests mean compactness values. For a xed number of superpixels, there are then much more ERS or SEEDS superpixels boundaries in the segmentation, which naturally increases the boundary recall value.

Undersegmentation error

Undersegmentation error (UE) is shown in Fig. 6(b). Given a ground truth and superpixel segmentation, this error measures the fraction of bleeding caused by superpixels that overlap a given ground truth segment. The standard formulation is

U E(s) = i k:s k ∩g i =∅ |s k -g i | i |g i | (8)
where s k are the outputs of the superpixel algorithm, g i the ground truth segments, and |.| denotes the size of an element.

There are signicant changes in this evaluation according to the authors, because it is not clear how to treat pixels that lie on a boundary between two labels. In [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF], the authors report a 5% tolerance margin, while authors in [START_REF] Liu | Entropy rate superpixel segmentation[END_REF] remove the boundaries of s k for the undersegmentation computation.

In this paper, we use the corrected undersegmentation error (CUE) proposed by the authors of SEEDS 6 dened as :

CU E(s) = k |s k -g max (s k )| i |g i | (9)
where g max (s k ) indicates the matching ground truth segment of s k with the largest overlap. This corrected undersegmentation error measure seems more accurate and does not relie on any ad-hoc solution.

Finally this fraction is simply averaged across all ground truth segments and all images. Fig. 6(b) shows that the undersegmentation error of ERGC is the lowest upon the considered stateoftheart methods.

Achievable Segmentation Accuracy

Achievable Segmentation Accuracy (ASA) is a performance upperbound measure. It gives the highest performance when taking superpixels as units for object segmentation. Each superpixel is labeled with the label of the ground truth segment with the largest overlap. The fraction of correctly labeled pixels is the achievable accuracy :

ASA(s) = k max i |s k ∩ g i | i g i (10)
6 http://arxiv.org/pdf/1309.3848v1.pdf

Compactness

We introduce this metric in addition to traditional ones to measure the compactness of the superpixels. Compactness represents the degree to which the superpixel shape is compact. It is dened as the ratio of the area of a region to the area of a circle with the same perimeter. It is calculated as follows:

COM P (s k ) = 4π|s k | p 2 k ( 11
)
where p k is the perimeter of the superpixel s k . The compactness is equal to 1 for a disc, π/4 for a square. 3), but produces more compact superpixels. Moreover, these performances are still higher than the popular algorithm SLIC.

Extensions

In this section, we propose several extensions to our initial framework, that illustrates the exibility of the approach.

3D extension

ERGC has been introduced for 2D images. It naturally extends to 3D volumes to produce supervoxels with minor modications of the algorithm.

Spatial constraint

For clarity of the display purposes, a spatial constraint has been added to F to generate images of Figure 3, 4 and 8. This constraint, similar to the one proposed in [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF], penalizes pixels p that are far from an initial seed s i , and is of the form p-s i 2 2

S

× m where m is the constraint parameter. In this case, the potential function is of the form:

F = C p -C i 2 2 + p -s i 2 2 S × m (12)
Adding such a spatial constraint increases the superpixels compactness and tends to produce square superpixels in at areas, see Figures 3,4. Nevertheless, as this is a constraint applied on the diusion, it decreases slightly BR, UE and ASA performances.

Combining ERGC with Edges

The proposed approach can easily be extended with additional terms to produce more powerful potential functions. We propose to combine the initial color-based potential function with edges maps. Based on structured forest, the approach proposed in [START_REF] Dollár | Structured forests for fast edge detection[END_REF] gives for each pixel p the probability E(p) that p belongs to an edge. Within the fast-marching algorithm, the potential of pixels belonging to region R i is computed as :

F = C p -C i 2 2 • + max γp i →s i E (13)
where is a small constant to avoid F being zero, and max γp i →s i (E) is the maximum edge probability along the geodesic γ p i →s i (see Figure 9).

Adding such an edge information to the potential is useful in case of microtextures such grass for instance. Edges probability is such areas is close to zero, making the potential of the pixels low. The propagating front is then not slowed down by the locally varying colors of the texture. Figure 11 compares the performances of this variant of ERGC, together with the initial ERGC method (without edges) and the top performer among the state-ofthe-art methods for each metric. Note that compacity values are not plotted for this method since there is no signicant changes from the initial ERGC method (Figure 6).

Iterative ERGC

The method proposed in this paper (detailed in Section 2) computes the superpixels in only one iteration. In this section, we go beyond this limitation and propose an iterative version of ERGC. Iterating the process allows to rene the seeds of the superpixels such that the next iteration produces better regions.

Given an initial superpixel R i with its seed s i , the renement is performed as follows :

Figure 9: Computation of the edge value of p i as the maximum edge probability along the geodesic γ p i →si .

1. Select all the pixels p i ∈ R i such that their color is the closest from C i , 2. from this set of seed candidates, select only those that are spatially the closest from s i . This renement scheme is illustrated in Figure 10. It may produce several new seeds per superpixels. In a case, these seeds share the same label, and form, after the diusion, a sole superpixel. The spatial selection is performed to avoid the generation of too many new seeds that may produce a degenerate superpixel after several iterations.

In our experiments, the number of iteration is xed to 10 and the conguration of seeds that lowers I U is retained.

At the cost of multiple iterations, this iterative variant enhances the whole performances of the algorithm (Figure 11).

For completeness purposes, Figure 11 also shows the performances of the iterative variant of ERGC with the use of edges maps. Finally, a comparison of relative processing times among all the methods is provided in Table 1. The pixels whose color is the closest from the superpixel color C i form a candidate set (gray dots, middle gure). From this set, only the spatially closest pixels (s 1 i and s 2 i ) are retained (right). Both s 1 i and s 2 i take the label of s i .

Conclusion

Superpixels have become an important preprocessing tool for many image based applications. In this paper, we proposed a method based on the eikonal equation that quickly creates accurate superpixels. Through the empirical experiments, we showed that ERGC outperforms existing superpixel methods in two of three metrics, while being outperformed by SEEDS and ERS on Boundary Recall. This last result has to be appreciated in the light of their superpixel compactness, the lowest within the state-of-the-art.

Proposed here for the generation of superpixels/supervoxels, we think that the diusion framework can oer an interesting choice in many other image processing tasks. ERGC can also be extended to arbitrary graphs to perform local or non-local image and data processing, and data clustering (supervertices).

Figure 1

 1 Figure 1 compares the geodesic distance map computed on a synthetic image (left) with the gradient-based potential function (middle) and the proposed one (right), with an initial seed depicted by the white dot. This example exhibits the main feature of the proposed potential function. With the gradient-based potential function (middle), the front propagates on the white square before having recovered all the black area. Some pixels belonging to the square then have a lower geodesic distance than pixels belonging to the black area. A good segmentation of the square based on these distances is then impossible. With the proposed potential function (right), the front propagates rst entirely on the black area before entering onto the white square. A good segmentation of the square can then easily be achieved.Such a result can be interpreted in the following way : For the gradientbased potential function, when a front arrives on a contour, its speed heavily

Figure 1 :

 1 Figure 1: Left : Initial image with the given seed depicted as the white dot. Middle : geodesic distance map U in false color obtained with the gradient-based potential function. Right : geodesic distance map obtained with the proposed potential function. A segmentation of the square based on the geodesic distances can not be obtained with the gradient-based potential function.

Figure 2

 2 Figure2compares the behavior of the proposed potential function and the gradientbased approach[START_REF] Desquesnes | Eikonal equation adaptation on weighted graphs: Fast geometric diusion process for local and non-local image and data processing[END_REF] on a color synthetic image with two seeds, and on a natural image with three seeds. In both cases, the proposed formulation gives a better segmentation.

Figure 2 :

 2 Figure 2: First row : initial image. The back and white circles depict the initial seeds. Second row : the result segmentation with the proposed F (left) and the gradientbased potential function (right). Third row : geodesic distances map U in fake color. Some isocontours are shown in white. Last row : segmentation of a natural scene. Seeds are depicted in black.

  2. add a new seed to the location of the maximum geodesic distance found at the previous step,[START_REF] Desquesnes | Eikonal equation adaptation on weighted graphs: Fast geometric diusion process for local and non-local image and data processing[END_REF]. let R i be the superpixel in which there is the new seed (red superpixel of Fig.3(b)), perform ERGC for pixels belonging to R i and its adjacent superpixels (blue superpixels of Fig.3(b). At this step, the seeds of the rened superpixels are left unchanged, 4. iterate steps 2 and 3 until the number of desired new seeds is reached.A resulting renement iteration, and the associated geodesic distance map are shown at Figures3(e) and 3(f ). The cost of a renement iteration depends on the number of pixels considered. Nevertheless, this renement is considerably less costly than a full ERGC pass, since it only deals with a small part of the image. Given superpixels of size S = N/K, and let be b the number of adjacent superpixels of R i , the complexity of one renement iteration is roughly O (bS log(bS)), which can be approximated with O N K log( N K ) . For information, one renement iteration costs approximatively 3 ms for 500 initial seeds, and 8 ms for 100 initial seeds on a Berkeley image. These time calculations have been obtained with a standard laptop equipped with an Intel mono core 1.30GHz processor and 4 GB RAM.

Figure 3 :

 3 Figure 3: (a) ERGC segmentation with 100 superpixels (approximatively). (b) Part of the image to be rened (red and blue superpixels). (c) Detail of the image (green box of (a)) and corresponding geodesic distance map (d). Adding a new seed to the location of the maximum distance results in a better segmentation (e). (f) New geodesic distance map corresonding to (e).

Figure 4 :

 4 Figure 4: Example of renement by moving superpixels. 3 superpixels have been removed

Figure 5 :

 5 Figure 5: Performances evolution with the renement steps : ERGC(a) stands for the initial pass of the algorithm, ERGC(b) after adding superpixels (Section 2.5), and ERGC(f) the nal results after adding and moving superpixels (Section 2.6). ERGC(m) shows performances of the algorithm with the spatial constraint m = 10 (equation 12).

Fig. 6 (

 6 Fig. 6(c) shows that ERGC gives the best Achievable Segmentation Accuracy compared to the other methods.

Fig. 6 (

 6 d) plots the mean compactness for each algorithm and for a dierent number of superpixels. Compactness for the grid of square superpixels (GRID) is the highest (around 0.85), and is not shown on Fig.6(d).

Figure 5

 5 Figure 5 plots the evaluation of the proposed algorithm and its evolution within the two proposed renement : the curves ERGC(a), ERGC(b)and ERGC(f ) stand for the algorithm performances without any renement, with only the adding superpixels renement, and with the two proposed renement respectively. One can particularly appreciate the improvements of the performances while rening iteratively the superpixels. Figure5also plots the performances of the proposed algorithm with a spatial constraint (m = 10 in equation 12) and SLIC (for comparison purposes). Adding a spatial constraint decreases slightly the BR, UE, and ASA performances of the algorithm in comparison to a potential based on color distances only (Equation3), but produces more compact superpixels. Moreover, these per-

Fig. 8

 8 Fig. 8 displays a supervoxels segmentation of a 3D volume from the MICCAI-2007 Grand Challenge dataset. For display purposes, only the interior part of the body is shown.

Figure 6 :

 6 Figure 6: Quantitative comparison of boundary recall (a), undersegmentation error (b), achievable segmentation accuracy (c), and compactness (d) with dierent numbers of superpixels.

Figure 7 :Figure 8 :

 78 Figure 7: Visual comparison of algorithms with 100 superpixels per image. First row: SLIC, second row: ERS, third row: SEEDS, fourth row: GrB, fth row: TurboPixels, sixth row: ERGC.

Figure 10 :

 10 Figure 10: Illustration of the iterative renement of the seed s i of a superpixel R i (left).

Figure 11 :

 11 Figure 11: Performances comparison between ERGC with the use of edges (ERGC(e)), iterative ERGC (ERGC(i)), iterative ERGC with edges (ERGC(ie)), ERGC, and the top performer among the state-of-the-art methods for each metric.

Table 1 :

 1 Relative processing times of several methods according to ERGC on a Berkeley dataset image, with 100 superpixels. These processing times may slightly vary with dierent numbers of superpixels. ERGC(e) and ERGC(ie) processing times include edge maps computation.

	Method	reference	Processing time factor
	SLIC	[1]	0.94
	ERS	[12]	5.2
	SEEDS	[18]	0.91
	GrB	[3]	0.95
	TP	[10]	20.2
	ERGC	This paper	1 (0.4s)
	ERGC(e) This paper	2.8
	ERGC(i) This paper	11
	ERGC(ie) This paper	13.8

Acknowledgements

This work is part of the SIRTDose project funded by the grant Physique Cancer 2012 (INSERM -Plan Cancer). The authors wants to thank the authors of [1] and [18] for the fruitlul discussions about superpixels and their evaluation. The authors also would like to thank the anonymous reviewers for their valuable comments and suggestions to improve the quality of the paper.