Ashley Hill
email: ashley.hill@cea.fr

Eric Lucet
email: eric.lucet@cea.fr

Roland Lenain
email: roland.lenain@inrae.fr

A new neural network feature importance method: Application to mobile robots controllers gain tuning

Keywords: Machine Learning, Neural Network, Robotics, Mobile Robot, Control Theory, Gain Tuning, Adaptive Control, Explainable Artificial Intelligence

This paper proposes a new approach for feature importance of neural networks and subsequently a methodology to determine useful sensor information in high performance controllers, using a trained neural network that predicts the quasi-optimal gain in real time. The neural network is trained using the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) algorithm, in order to lower a given objective function. The important sensor information for robotic control are determined using the described methodology. Then a proposed improvement to the tested control law is given, and compared with the neural network's gain prediction method for real time gain tuning. As a results, crucial information about the importance of a given sensory information for robotic control is determined, and shown to improve the performance of existing controllers.

Introduction

In robotic control, the search for more accurate and more adaptive controllers have been the foundation of research in the field. Where the controller must be able to adapt to varying known conditions or to react with respect to observed changes in the environment. However, it is not trivial to know a priori which aspects of the perception need to be included into the control law. As such many paths have been and are being explored for improving control law, such as computer vision [START_REF] Ha | World models[END_REF] or state estimators and observes [START_REF] Lenain | Robust sideslip angles observer for accurate off-road path tracking control[END_REF].

More recent papers have shown the use of neural networks for gain prediction [START_REF] Hill | Neuroevolution with cma-es for real-time gain tuning of a car-like robot controller[END_REF]. Where the gain changes with respect to the perception quality, allowing the controller to adapt to information that was underused. Indeed, due to their nature as universal function approximators [START_REF] Hornik | Universal approximation of an unknown mapping and its derivatives using multilayer feedforward networks[END_REF], neural networks can use most of the available sensor information; in order to predict a quasi-optimal gain that adapts to the changes in the perception.

However, many tasks cannot use neural networks, for safety reasons or in some cases for performance reasons. Indeed, they are considered black-boxes due to their mathematical complexity and high number of internal parameters [START_REF] Lecun | Deep learning[END_REF], making them hard to analyze and predict consistent behavior.

As such, the natural question that follows, is how to analyze trained neural controllers in order to understand their behavior, and possibly improve current control laws. This is an important question in the field of machine learning and artificial intelligence in general, and currently heavily researched [START_REF] Gunning | Explainable artificial intelligence (xai)[END_REF]. In this paper, the application of existing and novel analysis methods are used, in order to understand the how a trained gain prediction method reacts to its input information, and how this knowledge can be used in order to improve classic control law.

At first, a few methods for neural network analysis will be described, including a novel analysis method. Then, an experimental setup, where the gain prediction method is described, along with a dynamic simulation of a robotic car-like bicycle model. Followed by experiments showing how to integrate the analysis information into improving the control law. And finally, a discussion of the results and in depth analysis before concluding.

Analysis Method

Preliminaries

The analysis is based on feature importance, it describes how important each input feature is in order to obtain a good prediction. This is usually used in the context of decision trees [START_REF] Liaw | Classification and regression by randomforest[END_REF], where each node on the tree has a score determining the quality of its split. Which in some cases is the Gini impurity [START_REF] Suthaharan | Decision tree learning[END_REF]. The feature importance is described as the input parameters that lead to a low Gini impurity for each node that use the input feature for its split. When sorted, these feature importance show which inputs where the most useful in order to obtain a good prediction.

Unfortunately, decision trees struggle to outmatch the performance of neural networks, due to neural networks strengths as dimensional reducers and being universal function approximators for non-linear functions [START_REF] Lecun | Deep learning[END_REF]. This means that in most cases neural networks must be used in order to obtain the desired performance.

The notion of feature importance is still available to neural networks, however they are not as clear as for the decision trees. The most known method is the Temporal Permutation method, described in [START_REF] Molnar | Interpretable machine learning[END_REF], and detailed in the following section.

Feature Importance Using Temporal Permutation

Neural network predicts from a vector inputs, a desired output vector. By varying the inputs and observing the change in the output, a correlation can be established between each input and amount of change for the output. This correlation shows if a given input were to change, by how much would the output change.

If the assumption that the neural network predicts a quasi-optimal output is given. Then the change between the original predicted output, and the predicted output when the input was altered, should give the influence each input has on the output. And this is turn gives describes how important each input is to predicting the quasi-optimal output.

However, the changes that are made to each input must be in such a way, that the input values must remain consistent and realistic. For this, the approach proposed in [START_REF] Molnar | Interpretable machine learning[END_REF], use a temporal permutation method. Where for a list of input vectors recorded over time and a given input to analyze, the given input will be shuffled across all the input vectors over time. This causes a given input to be randomly permuted over time, meaning the input value will not longer hold any meaning with respect to the other inputs in the input vector. This input will then be similar to a random distribution of the same type as observed in the initial unshuffled input vectors.

When this method is applied to each input, a difference of the output, relative to each input can be achieved. Where the difference of the output is directly translated to the error of the output, if the assumption that the neural network predicts a quasioptimal output is given.

Feature Importance Using The Gradient

Feed forward multilayer perceptron neural networks, consist of a sequence of matrix multiplications, adds, and activation functions, from the given input to the given output [START_REF] Lecun | Deep learning[END_REF]:

y = a(b (n) + w (n,n-1) a(...b (1) + w (1,0) X))
Where y is the output, X is the input, a is the activation function, b (n) is the bias at the layer n, and w (n,n-1) is the weight matrix between the layer n and the layer n -1.

From this, the gradient between the output, and any component of the neural network can be achieved using the chain rule. Indeed this is the exact method that is used in backpropagation [START_REF] Lecun | Deep learning[END_REF] for gradient descent in supervised learning methods applied to neural networks.

However, backpropagation is only used to tune the parameter of the neural network in order to minimize an error between the predicted output and a desired output. A different method can be used with the chain rule to calculate the rate of change of the output value with respect to the input vector:

∂y ∂X = ∂y ∂a ∂a ∂z (n) ∂z (n) ∂s (n-1) . . . ∂s (1) ∂a ∂a ∂z (1) ∂z (1) ∂X ∂y ∂X = a (z (n) )w (n,n-1) . . . a (z (1) )w (1,0)
A variant of this method was previously used in image modification with neural networks, in order to change an input image that maximized a cost function [START_REF] Mordvintsev | Deepdream-a code example for visualizing neural networks[END_REF]. And for determining a Saliency map in image classifiers [START_REF] Simonyan | Deep inside convolutional networks: Visualising image classification models and saliency maps[END_REF].

Using this methods a jacobian matrix between each output component and each input component, can be obtained at each feed forward prediction of the neural network. With this, the average and variance of the rate of change for each input with respect to the output can be achieved over a given task. The contribution that allow the method to return the feature importance, is the assumption that the neural network predicts a quasi-optimal output. Where the rate of change of the output is directly translated to the error of the output. Meaning if a low rate of change for a given input is obtained, then the input does not contribute much to the quasi-optimal output, and as such is not considered to be an important feature for the prediction method.

Experimental Setup

Before any analysis of the neural network can be done, a training environment must first be established. For this, a robotic car-like model in a dynamic simulation is used to generate training samples. Which are then used as input to a covariance matrix adaptation evolution strategy (CMA-ES) method [START_REF] Hansen | The CMA evolution strategy: A tutorial[END_REF] in order to optimize a neural network's weights and biases, to minimize an objective function. This neural network has as a goal to output in real time, the quasioptimal gains for the steering controller, in a similar fashion to previous works [START_REF] Hill | Neuroevolution with cma-es for real-time gain tuning of a car-like robot controller[END_REF].

Robotic Model & Simulation

The robotic model is a dynamic bicycle model, which takes into account the slide slip angle and lateral forces applied to the rear and front axle. This model is show in the figure 1. The notation of the model being defined as follows: (D) is the trajectory being followed, s is the curvilinear abscissa along D, L is the wheel base length of the robot, v is the speed vector of the robot at the middle of the rear axle, θ is the angular error, y is the lateral error, c(s) is the curvature at the point s, δ F is the front steering angle, G is the center of mass of the robot, L R and L F are the distance from the center of mass to the rear and front axle respectively, F R and F F are the lateral force on the rear and front axle respectively, β is the vehicle sliding angle, α F and α R are the front and rear axle sliding angle respectively, I z is the moment of inertia across the Z axis.

(D) G θ δ F L L R L F α R α F β F R F F y v 1 c(s)
From this modeling, the following system of equations can be derived:

                             ṡ = v cos( θ) 1-c(s) y ẏ = v sin( θ) θ = 1 I z (-L F F F cos(δ F ) + L R F R ) β = -1 v 2 m (F F cos(β -δ F ) + F R cos(β)) - θ β R = arctan(tan β -L R θ v 2 cos(β) ) β F = arctan(tan β + L F θ v 2 cos(β) ) -δ F v 2 = v cos(β R ) cos(β)
The lateral forces applied to the rear and front axle, follow the Pacejka magic formula [START_REF] Bakker | Tyre modelling for use in vehicle dynamics studies[END_REF], in order to obtain a more realistic and dynamic environment.

An extended Kalman filter [START_REF] Welch | An introduction to the kalman filter[END_REF] is used in order to determine the robotic state from the sensor input, along with the covariance.

The steering actuators are modeled using an action delay of approximately 0.5s. This robot's steering is controlled using the following control equation:

δ F = arctan L cos 3 ε θ α k θ (e θ ) + κ cos 2 (ε θ ) with e θ = tan ε θ - k l ε l α
is the relative orientation error of the robot to reach its trajectory (i.e. ensuring the convergence of ε l to 0). This control detailed in [START_REF] Lenain | Robust sideslip angles observer for accurate off-road path tracking control[END_REF] guaranties the stabilization of the robot to its reference trajectory, providing a relevant choice for the gains k θ and k l .

Gain Prediction Model & Training

The gain prediction model used, is based on previous works [START_REF] Hill | Neuroevolution with cma-es for real-time gain tuning of a car-like robot controller[END_REF], where a neural networks predicts the quasi-optimal gain that minimize a given objective function. This neural network predict this gain using information that is underused in the control loop, such as the perception quality. This neural network is then optimized in order to minimize the objective function using the CMA-ES method [START_REF] Hansen | The CMA evolution strategy: A tutorial[END_REF]. The full control loop is shown on the figure 2: The objective function used is the following:

ob 1 = 1 T N ∑ n=0 |y(t n )| + L| θ(t n )| + k steer L|δ F (t n )| dt [m]
Where T is the total time taken to follow the path, N is the number of measured timesteps, t n is the time at the timestep n, y is the lateral error, θ is the angular error, δ F is the front steering, and dt is the time step between two samples. k steer is set to 0.5, as it showed to be the ideal compromise between minimizing the control errors, and keeping a low steering energy which minimizes oscillations.

Feature Importance for robotic control

In robotics, there can be many sensors that can measure a wide variety of different kinds information. However, it is not always clear what kind of information or sensors will be of use in order to develop a performant control equation.

As such, the goal of the following experiments, is to use a trained neural network that can predict the quasi-optimal gain. And from this, determine which inputs of the neural network are used the most in order to minimize the given objective function.

This will give a list of features that are needed in order to improve the performance of the control equation. Some of which will then be integrated into the control law.

Results

The following results were obtained using the previously described methods, with a line following task and the trajectory shown in the figure 3, at 2.0m.s -1 , 1.5m.s -1 , and 1.0m.s -1 . Midway though the trajectory, a GPS noise of 1m is applied to simulate a perception quality loss. Two change lanes occurs along the trajectory, in order to simulate a GPS constellation jump or a sudden change in the target set-point.

Baseline

The baseline method that will be used to compare in the following sections, is the expert tuned constant gain method. Using a gain of k l = 0.2, and k θ = 1.0.

The results for the experiment can be see in the figure 4. Where the baseline method reached an end objective function value of 0.389m, 0.336m, and 0.308m for 2.0m.s -1 , 1.5m.s -1 , and 1.0m.s -1 respectively. This method has some obvious shortcomings, as it is not adaptive to the changes in speed or sensor accuracy. Which can be observed in the noisy region, where the controller is reacting to the GPS noise.

Gain Adaptation

The neural network gain method uses the information in the control loop, in order to minimize the objective function. As such, it is able to predict the quasioptimal gain along the trajectory at any given time.

The results for the experiment can be see in the figure 4. Where the neural network gain method reached an end objective function value of 0.372m, 0.325m, and 0.291m for 2.0m.s -1 , 1.5m.s -1 , and 1.0m.s -1 respectively. The neural network is adapting the gain with respect to changes in the speed, sensor accuracy, curvature, and error. This allows the method to lower it's objective function substantially when compared to the baseline method.

The feature importance can now be done over the trained neural network. For this, both the temporal permutation feature importance and the novel gradient feature importance are done.

On figure 6, the temporal permutation feature importance can be observed. It is shown as the absolute mean difference between the original predicted gain, and the predicted gain when the given input is shuffled over time. From this, the most important features from highest to lowest are the Kalman covariance matrix denoted C xy , the speed and target speed denoted v and v target respectively, the lateral error denoted y, the angular error denoted θ, and the curvature denoted c(s).

The temporal permutation feature importance describes the absolute change in the gain, with respect to the input feature. However is it does not show the rate of change of the gain, with respect to the values of the input features. For this the gradient feature importance is used.

On figure 7, the gradient feature importance can be observed. It is shown as the mean and the variance of the jacobian matrix of the gain for each input. Sim- ] ilarly to the temporal permutation feature importance, the most important features from highest to lowest are the Kalman covariance matrix denoted C xy , the speed and target speed denoted v and v target respectively, the lateral error denoted y, the angular error denoted θ, the curvature denoted c(s), and the future curvature (the predicted curvature at t + 1s) denoted future c(s) due to its high variance.

Cxy [m 2 ] vtarget [m. s 1 ] dy/dt [m. s 1 ] /dt [rad. s 1 ] future c(s) [m 1 ] f [rad] r [rad] v [m. s 1 ] 2.50 2.25 0.5 0.0 0.5 k -mean jaccob matrix of dk dx y [m] [rad] c(s) [m 1 ] Cxy [m 2 ] vtarget [m. s 1 ] dy/dt [m. s 1 ] /dt [rad. s 1 ] future c(s) [m 1 ] f [rad] r [rad] v [m. s 1 ] 4.50 4.25
From this, the equivalence of the feature importance methods can be implied. However the gradient feature importance has some strong strengths to it: It does not need each input features to vary in order to obtain the feature importance. Indeed if a given input does not explore the span of values during the analysis, the temporal permutation feature importance will not return the correct feature importance. Furthermore the gradient feature importance can return the jacobian matrix for each input and output, at each feed forward of the neural network. Allowing for real time analysis and approximate prediction of the behavior of the neural network. And finally, the mean jacobian matrix for each input and output can be used, as an approximation of the neural network's behavior, that can be exploited to improve the control equation, as shown in the following section.

Improving Control Law

The ideal way to change the control law, is to derive it from the model while taking some of the important features into account (covariance, speed, sliding angles, ...). However, in this case study a simple modification of the gain will be used in order to show a proof of concept.

For this, the gains equation are augmented using the mean jacobian ∂y ∂X matrix for the speed. In this case study the covariance is not used, due to how the neural network is using it for non-linear adaptation of the gain. As such a linearization of the jacobian is not a valid approximation of the original gain behavior for the covariance, and in which during experimentation lead to unstable behavior.

Using the mean jacobian matrix, the following gain equations are derived:

k l = v ∂ kl ∂v + k lv k θ = v ∂ kθ ∂v + k θv
Where ∂ kl ∂v , ∂ kθ ∂v are the mean jacobian of the predicted gains with respect to the speed, v is the longitudinal speed, and k lv , k θv are the new tune-able gains for the control equation. This modification should allow the control law to change its reactivity to the error with respect to changes in the speed.

The results for the experiment can be see in the figure 8. Where the improved control equation method reached an end objective function value of 0.387m, 0.327m, and 0.288m for 2.0m.s -1 , 1.5m.s -1 , and 1.0m.s -1 respectively. The improved control is adapting the gain with respect to changes in the speed. This allows the method to lower it's objective function substantially when compared to the baseline method, and to achieve similar performance to the neural network gain adaptation method. Furthermore, it is able to capture the fast initial convergence to the trajectory that the neural network gain adaptation method had, thanks to the initial speed up (visible from t = 0 to t = 5).

1.0m.s -1 1.5m.s -1 2. In the table 1, the objective function for each method and speed can be observed. In all cases the baseline constant gain method had the highest objective function, meaning it had the worst performance. For the improved control and the neural network gain method, both results are comparable, as most of the performance gained is thanks to the speed adaption. However the improved control does not adapt to changes in the covariance, which in some cases allows the neural network gain method to outperform the improved control method.

A novel method for feature importance and a novel methodology to determine useful sensor information was proposed. This feature importance method allows the analysis of a neural network's behavior, to show the importance of each sensor information, and to potentially build a linear approximation of the neural network for a given input.

It has been applied to a steering controller of a carlike robot for a line following task in a highly dynamic simulated environment. In order to analyze a gain prediction method, and determine the optimal changes to the control equations to improve its performance. Indeed, the tested modification to the control law has been shown to reach comparable performance to the initial neural network gain prediction method.

This methodology can be applied to any given simulated model of a robot control task, in order to improve its control performance for a given criteria encoded as a objective function.

However, the sensor information must be ideally used to derive new control law from the robotic model, as using a linear approximation for a neural network will not encode the complete characteristic behavior to the neural network. As such this methodology far more powerful as a tool to describe what is important for control law, not how to derive a novel control law.

Future works include validating this methodology on varying control tasks in different field, and to use the novel feature importance method to assist in demystifying neural networks.

Figure 1 :

 1 Figure 1: The dynamic robot model.

  Figure 3: The trajectory

Figure 4 :

 4 Figure 4: Above: The method line following. Below: the objective function over time.

Figure 5 :

 5 Figure 5: Above: The method line following. Middle: The predicted gain, where the reference gain is the baseline constant gain. Below: the objective function over time.

Figure 6 :

 6 Figure 6: The temporal permutation feature importance. Above: for the k l gain. Below: for the k θ gain.

Figure 7 :

 7 Figure 7: The gradient feature importance, as the mean jacobian matrix. The figure is cut midway to avoid scaling issue of outliers. The error bars, show the variance of the jacobian matrix. Above: for the k l gain. Below: for the k θ gain.

Figure 8 :

 8 Figure 8: Above: The method line following. Middle: The predicted gain, where the reference gain is the baseline constant gain. Below: the objective function over time.

Table 1 :

 1 The objective function obtained for each method and speed over the trajectory.

	0m.s -1


Eric Lucet
email: eric.lucet@cea.fr

Alain Micaelli
email: alain.micaelli@cea.fr

A kinematic parameter estimator applied to a bi-directional vehicle

Keywords:

This paper deals with the online estimation of the geometric and kinematic parameters of a wheeled mobile robot, with the objective of its precise navigation. To do this, the implementation of an algorithm for estimating these parameters is proposed. It is based on a quadratic criterion to be minimized, which is a function of the difference between the measured and estimated robot pose. The estimator inputs external measurements of robot positions and speeds and proprioceptive inertial and odometer measurements, and outputs an estimate of model parameters. Its expression does not depend on the position of the robot relative to the path to be followed and its possible tracking errors, nor on the speed of the tracking. The experimental implementation carried out on a real bi-directional container truck under realistic operating conditions has shown its performance.

INTRODUCTION

One of the main challenges of mobile wheeled robots is their precise navigation along a trajectory. It depends among other things on their type of wheels and their configurations (see [START_REF] Tzafestas | Introduction to Mobile Robot Control[END_REF] for an overview of wheels and chassis). This has an impact on the models of the robot to be considered and their parameters. Possible models are, for example, sensor models, actuator models, kinematic models, dynamic models or flexibility models. And an accurate identification of a robot model from measurements is very useful, either for state observation [START_REF] Carvalho Filho | The Impact of Parametric Uncertainties on Mobile Robots Velocities and Pose Estimation[END_REF]) or for the implementation of a model-based control of the robot (see Fig. 1). State observers can use Bayesian estimation theory tools such as the Extended Kalman Filter (EKF) or the Particle Filter [START_REF] Thrun | Probabilistic Robotics[END_REF] or [START_REF] Siegwart | Introduction to Autonomous Mobile Robots[END_REF]), while many model-based mobile robot controllers [START_REF] Frank | Control Theory Tutorial: Basic Concepts Illustrated by Software Examples[END_REF]) are available, such as model reference adaptive control [START_REF] Ashoorirad | Model reference adaptive path following for wheeled mobile robots[END_REF]) or model predictive control (MPC) [START_REF] Rossiter | Model-based predictive control: a practical approach[END_REF], [START_REF] Lenain | Adaptive and predictive non linear control for sliding vehicle guidance: application to trajectory tracking of farm vehicles relying on a single rtk gps[END_REF]).

The state reconstructor or estimator is a system whose inputs are the inputs and outputs of the real process and whose output is an estimate of the state of that process. The process model can be linear or non-linear, continuous or discrete, deterministic or stochastic. The case of noisy systems, where random phenomena occur, is referred to filters such as EKF. For example, we can cite the case of an EKF-based kinematic parameter estimation for a passenger vehicle [START_REF] Brunker | GNSS-shortages-resistant and self-adaptive rear axle kinematic parameter estimator (SA-RAKPE)[END_REF]). Other deterministic methods have been proposed in the past. Non-linear model estimation can be done with the Gauss-Newton method [START_REF] Norton | An Introduction to Identification[END_REF]) which ensures quadratic convergence. In that case, the estimation is done iteratively by minimizing the sum of the squares of errors between the measurement outputs and the outputs predicted by the current theoretical model. Since input mea-surement noise can lead to bias errors, this constraint can be processed with total least squares [START_REF] Van Huffel | The Total Least Squares Problem: Computational Aspects and Analysis[END_REF]). Less accurate, the screw axis measurement [START_REF] Hollerbach | The Calibration Index and Taxonomy for Robot Kinematic Calibration Methods[END_REF]) is a class of methods which measure the joint axes in the form of lines in space. One of these methods is the circle point analysis [START_REF] Mooring | Fundamentals of manipulator calibration[END_REF]) for the moving of joints one by one to generate a circle around a measurement point. Other approaches use the Jacobian matrix [START_REF] Bennett | Kinematic calibration by direct estimation of the Jacobian matrix[END_REF]), or even rely on the control theory. One example is the estimation of slip parameters using the sliding-mode control theory [START_REF] Song | Non-linear observer for slip parameter estimation of unmanned wheeled robots[END_REF]), based on the vehicle kinematic model and measurements.

Here, the deterministic case is investigated for the observation of the parameters of a discretized non-linear model of a mobile robot, for its implementation on the embedded computer of a real robotic system. The proposed new approach is similar to those mentioned above in that a quadratic criterion is minimized. Nevertheless, it has the advantage of having an expression that does not depend on trajectory tracking performance, and its expression is sufficiently general to be able to propose the estimation of any type of geometric or kinematic parameter of any type of wheeled robot model. The block diagram Fig. 1 gives an overview of the classic architecture of the control system (similar to the one of [START_REF] Cariou | Automatic guidance of a four-wheel-steering mobile robot for accurate field operations[END_REF]).

State observer

Controller

Mobile Robot Trajectory

Model parameters Estimator

Fig. In particular, this parameter estimation makes it possible to compensate for wheel-ground sliding phenomena with dynamics dependent on the position measurement period, but also to estimate the value of parameters that are difficult to measure directly with precision, such as the effective wheel diameter (required for longitudinal odometry), steering angle offsets, or even the wheelbase.

This architecture ensuring an online correction of the model parameters, this will make it easier to manage high speed variations.

Thereafter, the kinematic model of a vehicle with steered front and rear axles is designed in section 2 for observation and control. The model chosen is horizontal plane -roll and pitch are not taken into account. Then, an estimator of the parameters of this model is designed in section 3. This estimator is applied to particular cases of parameters in section 4. Then, test results are presented in section 5 before the conclusion section 6.

KINEMATIC MODEL

The kinematic model used for the estimator is described in this section. The plan model of a bi-steered vehicle is chosen, it includes the usual case of the car. However, any other type of model could have been considered, the method remaining the same.

Equations of motion

The kinematic model is computed by using the Lie group formalism, described in detail for example by [START_REF] Murray | A Mathematical Introduction to Robotic Manipulation[END_REF]. This allows more clarity in equations writing.

In the following, The frame (f ) is fixed to the vehicle front wheel centre with its longitudinal axis in the wheel steering direction. The same applies to frame (r). ( 0) is the world frame and (m) is the vehicle body frame. (c f ) and (c r ) are the front and rear reference frames on the path. Notations are the following. j T k i is the twist of frame (i) with respect to frame (j), reduced at the origin of frame (k). For simplicity, (k) is not specified if k = i. Notations are the same for the pose matrix j H i . j Ad i is the adjoint matrix, for transporting frame (i) into frame (j). Notations are the same for j R i , j p i and scalars.

In the Fig. 2, α f,r are the front and rear steering angles, δ f,r are the front and rear sliding angles and β f,r = α f,r + δ f,r . Then, kinematic equations are defined in both frames (f ) and (r) by the following twists:

o T f =   θm + βf v f 0   , o T r =   θm + βr v r 0   (1)
And twists of frames (f ) and (r) relative to the vehicle body and expressed in these frames are written as follows:

Fig. 2. Horizontal plane model of a bi-directional vehicle relative to its reference trajectory

m T f =   βf 0 0   , m T r =   βr 0 0   (2)
The twists defined by ( 1) are not independent. The combination of ( 1) and ( 2) allows to link the longitudinal velocities v f and v r in order to obtain the following expression of the yaw velocity θm of the mobile:

o T m = m Ad f ( o T f -m T f ) = m Ad r ( o T r -m T r ) (3)
with:

• m Ad f,r = 1 0 1×2 -m p ⊥f,r m R f,r ;
• m p f,r , the vectors mf and mr expressed in the mobile frame (m); • m p ⊥f,r , the vector obtained by a π 2 rotation of the vector m p f,r ; • m R f,r , a rotation matrix defining the direction of frames respectively (f ) and (r) in the (m) frame.

From the above equations, the following two equations are obtained:

θm = v f σ f v r = v f φ v f (4)
with following variables:

• l = l f -l r (l f and l r are signed with respect to the (m) frame)

• σ f = sin(β f -βr) l cos βr • φ v f = cos β f cos βr
Similarly, taking v r , βf and βr as independent variables leads to:

θm = v r σ r v f = v r φ v r (5)
with following variables:

• σ r = sin(β f -βr) l cos β f • φ v r = cos βr cos β f

Estimation models

For the estimation, twists and positions of (f ) and (r) frames are used. Differently from previous works (e.g. [START_REF] Cariou | Automatic guidance of a four-wheel-steering mobile robot for accurate field operations[END_REF], [START_REF] Lenain | Adaptive and predictive non linear control for sliding vehicle guidance: application to trajectory tracking of farm vehicles relying on a single rtk gps[END_REF] or [START_REF] Lenain | Sideslip angles observer for vehicle guidance in sliding conditions: application to agricultural path tracking tasks[END_REF]), the parameters estimation does not depend on a path. The major advantage is that it is not disturbed in case of significant path tracking error.

MODEL PARAMETERS ESTIMATOR ALGORITHM

The objective of this section is to define an estimator of geometric and kinematic parameters, based on the kinematic model defined in the previous section and by defining the variation of the parameter to be estimated such that it minimizes a quadratic criterion that is a function of the difference between the measured and estimated robot pose.

For the estimator design, measured proprioceptive state variables are considered, such as:

• ω f and ω r , the averages of speeds of rotation of the front and rear wheel axles; • α f and α r , the front and rear steering angles;

• αf and αr , the front and rear steering velocities.

These odometry measurements are received from the vehicle at the sampling period dt. Inertial data could also have been considered. In practice, the controller's steering speed outputs are often used because they avoid measurement noise.

For the exteroceptive variables measurements, the following variables are assumed to be available: positions of (f ) and (r) frames, are given by the following relationships:

• o p mmes ,
         o Tf,r = m Ad -1 f,r o Tm + m Tf,r o Ĥf,r = o Ĥm +   βf,r l f,
where o Tk+j k+j odo is the odometry twist calculated at the k + j step. It is important to note that, in general, the product of matrix exponentials is not equal to the exponential of the sum of these matrices, and in particular: ) .

The expression of o

Ĥf,r k+n can also be written as follows:

o Ĥf,r k+n = o R f,r k R (ω k dt) • • • R (ω k+j dt) • • • o p f,r k + o R f,r k f,r k p f,r k+1 + • • • .
Based on a variation of a parameter µ to identify, the variation of o Ĥf,r k+n expressed in frames (f ) and (r) by the equation ( 7), may be written as:

∂ ∂µ o Ĥf,r k+n = ∂ ∂µ o Ĥf,
A criterion C to be minimized is then considered:

Min µ C = 1 2 o Ĥf,r k+n -1 [ o H f,rmes ] -o Ĥf,r k+n v 2
To simplify the notations, it is rewritten as follows:

Min µ C = 1 2 [H] -1 δ [H] v 2
And its derivative is calculated:

∂ ∂µ [H] -1 ([H mes ] -[H]) = ∂ ∂µ [H] -1 [H mes ] = -[H] -1 ∂ [H] ∂µ [H] -1 [H mes ] = -[H] -1 ∂ [H] ∂µ Id + [H] -1 δ [H]
For a small variation of µ:

∂ ∂µ [H] -1 δ [H] δµ = -[H] -1 ∂ [H] ∂µ δµ -[H] -1 ∂ [H] ∂µ [H] -1 δ [H] δµ
The second term of the right member is considered null for a second order variation. Therefore, the variation of criterion C for a variation of µ is:

δ µ C = [H] -1 δ [H] vt ∂ ∂µ [H] -1 δ [H] v δµ ≈ -[H] -1 δ [H] vt [H] -1 ∂ [H] ∂µ v δµ
To guarantee a negative or null value of this derivative, δµ is chosen to be:

δµ = +l µ [H] -1 ∂ [H] ∂µ vt [H] -1 δ [H] v = +l µ [H] -1 ∂ [H] ∂µ vt f,r k+n Tf,r k+n f,r k+n mes v
where l µ is a positive factor. This last equality is clarified:

[H] -1 δ [H] = R t -R t p 0 1x2 1 δR δp 0 1x2 0 At first order: R t δR = R t R mes -Id 2x2 = cos (δθ) -sin (δθ) sin (δθ) cos (δθ) -Id 2x2 ≈ δθJ with J = 0 -1 1 0 . Then: [H] -1 δ [H] ≈ δθJ R t δp 0 1x2 0 = [T δ ]
with T δ = δθ R t δp being f,r k+n Tf,r k+n f,r k+n mes here.

Finally, by using equation ( 8), the following result is obtained:

   δµ = l µ dt 1≤i≤n ∂ ∂µ o Tk+i t k+i odo k+i Ad -t k+n f,r k+n Tf,r k+n f,r k+n mes μ+ = μ-+ δµ (9)
The estimator ( 9) is applied to some µ parameter cases, namely wheel diameter, wheelbase and sideslip angle, in the following section.

PARAMETER ESTIMATES

Wheel diameter estimate

Estimating the equivalent wheel diameter is equivalent to estimating the longitudinal slippage. For estimating the equivalent wheel diameter, equations ( 1) are considered.

To simplify them, βf and βr are supposed to be invariant during a sampling step, so odometry twists o Tf,r j are simplified accordingly.

Then, by using equations ( 4) and ( 5), and that v f,rj = d wheel 2 ω f,rj , the partial derivative of the odometry twist o Tf,rj with respect to d wheel is given by:

∂ o Tf,rj ∂d wheel = 1 2   ωf,rj σf,rj ωf,rj 0   (10)
Depending on whether the front or rear wheels of the mobile robot are actuated, the expression with o Tfj or o Trj will be used. The Estimator (9) can then be used with a gain l d wheel .

Wheelbase estimate

For estimating the mobile robot wheelbase, odometry twists o

Tf,r j are simplified as previously, with βf and βr being supposed to be invariant during a sampling step. And, by using equations ( 4) and ( 5), the partial derivative of the odometry twist o Tf,rj with respect to l is given by:

∂ o Tf,rj ∂l =    - 1 l vf,rj σf,rj 0 0    (11)
The Estimator ( 9) is then used with a gain l l .

Sideslip angles estimate

For estimating the sideslip angles, βf and βr are also supposed to be invariant during a sampling step, and sideslip angles are now taken into account. Thus, odometry twists are rewritten differently from equation ( 1) as follows:

o Tf,rj =    θm vf,rj cos δf,rj vf,rj sin δf,rj    (12)
Considering an actuation of the front wheels by using equations (4), partial derivatives of odometry twists o Tf,r j with respect to sliding angles δ f,r are given by:

                 ∂ o Tfj ∂δ f =   vfj ψf fj 0 vfj   ∂ o Trj ∂δ r =    vfj -ψf fj + φr fj σfj 0 vfj φv fj    (13)
with:

• ψ f f = cos(β f -βr) l cos βr ; • φ r f = tan β r ;
Considering an actuation of the rear wheels by using equations ( 5), partial derivatives of odometry twists o Tf,r j with respect to sliding angles δ f,r are given by:

                 ∂ o Tfj ∂δ f =    vrj -ψr rj + φf rj σrj 0 vrj φv rj    ∂ o Trj ∂δ r =   vrj ψr rj 0 vrj   (14)
with:

• ψ r r = -cos(β f -βr) l cos β f ; • φ f r = tan β f ;
Terms cos δ f,r j and sin δ f,r j are simplified in equations ( 13) and ( 14).

Estimator ( 9) is then used with the gains l δ f and l δr for the estimation of respectively δ f and δ r .

CALIBRATION OF A BI-STEERABLE CONTAINER TRUCK

Description of the test

If the parameter estimator can be used to estimate online parameters during vehicle navigation, it can also be used to estimate steering offsets. This is what has been done for the front and rear directions of an electric prototype container truck, designed for commercial port operations.

It is an 18tons empty truck, with a 2m track and an 11m wheelbase. The actuation of the front and rear steering axles is hydraulic, by flow control servovalves. The lowlevel control is a simple loop proportional to the steering angle error. Since the front wheels are electrically actuated, equations (13) are used.

Proprioceptive sensor data are angular increments of the sensors on both front wheels and the crown steering angle of front and rear axles. Those data are received from the vehicle low-level controller through a CAN bus every 10ms. Exteroceptive sensors used for this test are a lidar and six pairs of ultrasonic distance sensors (24ms response time, 0.069mm resolution, range between 20mm and 150mm). The measured position and orientation of the truck are obtained from the data of these sensors after filtering. This test requires a high degree of accuracy and is performed inside a docking station where containers are placed and deposited by trucks. Outside this docking station, the laser allows the truck to be positioned relative to it. Inside, ultrasound is used to position the truck in relation to the flat walls of this station. The truck moving at low speed, and the ultrasonic sensors being already almost orthogonal to the wall, this allows an accurate millimetric measurement. Then, the steering sensor calibration process is performed using the sideslip angle estimator while the truck moves in a straight line for 30s with zero front and rear steering angle setpoints, at a low speed of 0.25m/s. Tests are carried out under realistic conditions on asphalt roads.

Results analysis

The position of the truck is shown in Fig. 3, with a solid black line for the position measured in the centre of the vehicle. The orientation measured at the centre of the vehicle is shown in Fig. 4 in solid black line. It is noticeable that the position is quite noisy with an amplitude of about 2 to 3cm. This is mainly due to the long length of the truck, which makes it difficult to position it accurately, as any change in the orientation of the vehicle has an impact on its lateral positioning. The walls of the docking station, which are not perfectly flat, can also deteriorate the quality of ultrasound measurements.

As for the measured orientation, although almost zero, it fluctuates slightly between 0.016rad (0.9deg) and 0.026rad (1.5deg) with a progressive increase until 20s before decreasing again. These fluctuations have a significant influence on the estimation of sideslip angles, with the measured position and orientation terms being used in the ex-pression of the twist f,r k+n Tf,r k+n f,r k+n mes used in equation ( 9).

The sideslip angles are shown in Fig. 5, as a cyan solid line for the front axle and a dotted blue line after filtering, an orange solid line for the rear axle and a red solid line after filtering.

Fig. 5. Sideslip angles

The front and rear sideslip angles take about 10s to converge. The value reached is about 0.01rad for both axles. Then, a variation that is maximum at t = 20s is observed, then again the curves converge towards the value of 0.01rad. This clearly shows the dependence of the algorithm on measurement noise and quality, and here in particular on the orientation measurement of the vehicle o θ mmes .

CONCLUSION

An estimator of the geometric and kinematic parameters of wheeled mobile robots was presented. This estimator was applied to the particular case of a bi-steered vehicle model, although it would also be appropriate for any other type of kinematics, such as the multibody system presented by [START_REF] Lucet | Stabilization of a roadtrain of articulated vehicles[END_REF]. It allows an online update of control or observation models, or a calibration of odometric sensors. This algorithm has the advantage of being independent of the trajectory to be followed, and it is insensitive to tracking errors.

Experimental results have demonstrated its effectiveness in accurately estimating parameters, with a reasonable convergence time of a few seconds. However, it remains sensitive to exteroceptive measuring noises of the position and orientation of the vehicle, which therefore remain to be verified for the correct use of this estimator.

An extension to parameter estimation of a dynamic model is envisaged, although additional difficulties of dependence on the speed of the robot have to be addressed.

  the measured position of frame (m) with respect to the reference frame (o); • o θ mmes , the measured mobile robot orientation with respect to the reference frame (o). The state observation is based on equation (3) and on the pose o H m = o θ m o p m of the mobile. Observed speeds and

  of the right member is null since µ is varying only after the step k. Then, by rearranging the first product:

Fig. 3 .

 3 Fig. 3. Truck path tracking

  1. Control block diagramThe Controller and the State observer are based on a kinematic model of the vehicle. The function of the Model parameters Estimator is to correct the parameters of this model online, from the controllers' control outputs and the state observer's position and speed outputs.


Stéphanie Jehan-Besson

Muriel Gastaud

Frédéric Precioso

Michel Barlaud
email: barlaud@i3s.unice.fr

Gilles Aubert

Éric Debreuve

From snakes to region-based active contours defined by region-dependent parameters

This synthetic paper deals with image and sequence segmentation when looking at the segmentation task from a criterion optimization point of view. Such a segmentation criterion involves so-called (boundary and region) descriptors which, in the general case, may depend on their respective boundary or region. This dependency must be taken into account when computing the criterion derivative with respect to the unknown object domain (defined by its boundary). If not, some correctional terms may be omitted. This article focuses computing the derivative of the segmentation criterion using a dynamic scheme. The presented scheme is general enough to provide a framework for a wide variety of applications in segmentation. It also provides a theoretical meaning to the active contour philosophy.

Introduction

This paper does the synthesis of several years of development in active contour segmentation conducted at I3S (Informatique, Signaux et Systèmes de Sophia Antipolis) and Dieudonné laboratories, CNRS (French National Research Center) and University Of Nice-Sophia Antipolis, France.

The purpose of segmentation is to isolate an object (or several objects) of interest in an image or a sequence. Given an initial contour (a closed curve), the active contour technique consists in applying locally a force (or displacement, or velocity) such that the initial contour evolves toward the contour of the object of interest. This force is derived from a characterization of the object formally written as a criterion to be optimized.

A. Boundary-Based Active Contours

In boundary-based active contour techniques, the object is characterized by properties of its contour only. The original active contour developments were called snakes. [START_REF] Kass | Snakes: Active contour models[END_REF] Only the convex hull of objects could be segmented because these techniques were based 1 on a minimum length penalty. In order to be able to segment concave objects, a balloon force was heuristically introduced. It was later theoretically justified as a minimum area constraint 2 balancing the minimum length constraint. The geodesic active contour technique [START_REF] Caselles | Geodesic active contours[END_REF] is the most general form of boundary-based techniques. The criterion corresponding to this technique is

J(Γ) = Γ k(x) dx ( 1
)
where Γ is a contour, k is a positive function "describing" the object of interest, and

x is a point of the image. The active contour evolution equation [START_REF] Caselles | Geodesic active contours[END_REF] is

∂Γ ∂τ = (κk -∇k • N )N ( 2
)
where τ is the evolution parameter (Γ(τ = 0) = Γ 0 , initial contour, and Γ(τ → ∞) → segmentation), κ is the curvature of Γ, operator • represents the inner product, and N is the inward normal. Actually, this equation should be written with κ(x), k(x), and N (x) where x is a point of Γ(τ ). The contour minimizing J can be interpreted as the curve of minimum length in the metric defined by function k. If k is the constant function equal to one, Eq. ( 2) is called the geometric heat equation by analogy with the heat diffusion equation. Function k can also be a function of the gradient of the image. For instance,

k(x) = 1 1 + |∇I(x)| (3)
where I is the image and ∇ is the intensity gradient. In this case, the object contour is simply characterized by a curve following high gradients. As a consequence, the technique is effective only if the contrast between the object and the background is high. Moreover, high gradients in an image may correspond to the boundaries of objects that are not of interest. Regardless of function k, information on the boundary is too local for segmentation of complex scenes. A global, more sophisticated object characterization is needed.

B. Region-Based Active Contours

In order to better characterize an object and to be less sensitive to noise, region-based active contour techniques were proposed. [START_REF] Cohen | Reconstruction of digital terrain model with a lake[END_REF][START_REF] Ronfard | Region-based strategies for active contour models[END_REF] A region is represented by mathematical expressions called "descriptors" in this paper. Two kinds of region are usually considered: The object of interest and the background. Note that region-based/boundarybased hybrid techniques are common. [START_REF] Chakraborty | Deformable boundary finding in medical images by integrating gradient and region information[END_REF][START_REF] Zhu | Region competition: unifying snakes, region growing, and bayes/MDL for multiband image segmentation[END_REF][START_REF] Paragios | Geodesic active regions for motion estimation and tracking[END_REF][START_REF] Paragios | Geodesic active regions and level set methods for supervised texture segmentation[END_REF] In the general case, descriptors may depend on their respective regions, for instance, statistical features such as mean intensity or variance within the region. [START_REF] Yezzi | A statistical approach to snakes for bimodal and trimodal imagery[END_REF] The general form of a criterion including both regionbased and boundary-based terms is

J(Ω in ) = Ω in k in (Ω in , x) dx + η Γ k b (x) dx (4)
where Ω in is the domain inside Γ and η is a positive constant. An example of descriptor

k in is k in (Ω in , x) = (µ(Ω in ) -I(x)) 2 (5)
where µ is the mean intensity

µ(Ω in ) = Ω in I(x) dx Ω in dx . ( 6
)
Function ( 3) is an example of descriptor k b . Classically the integral on domain Ω in in criterion ( 4) is reduced to an integral along Γ using the Green-Riemann theorem [START_REF] Chakraborty | Deformable boundary finding in medical images by integrating gradient and region information[END_REF][START_REF] Zhu | Region competition: unifying snakes, region growing, and bayes/MDL for multiband image segmentation[END_REF][START_REF] Paragios | Geodesic active regions for motion estimation and tracking[END_REF][START_REF] Chesnaud | Statistical region snake-based segmentation adapted to different physical noise models[END_REF][START_REF] Samson | A level set model for image classification[END_REF][START_REF] Chan | Active contours without edges[END_REF] or continuous media mechanics techniques. [START_REF] Debreuve | Space time segmentation using level set active contours applied to myocardial gated SPECT[END_REF][START_REF] Amadieu | Inward and outward curve evolution using level set method[END_REF] Two active contour approaches are possible to minimize the resulting criterion: (i) It is possible to determine the evolution of an active contour from τ to τ + dτ without computing a velocity: The displacement of a point of the contour is chosen among small random displacements as the one leading to the locally optimal criterion value. [START_REF] Chakraborty | Deformable boundary finding in medical images by integrating gradient and region information[END_REF][START_REF] Chesnaud | Statistical region snake-based segmentation adapted to different physical noise models[END_REF] However this implies to compute the criterion value several times for each point; (ii) Alternatively, differentiating the criterion with respect to τ allows to find an expression of the appropriate displacement, or velocity, for each point. [START_REF] Zhu | Region competition: unifying snakes, region growing, and bayes/MDL for multiband image segmentation[END_REF][START_REF] Paragios | Geodesic active regions for motion estimation and tracking[END_REF][START_REF] Samson | A level set model for image classification[END_REF][START_REF] Chan | Active contours without edges[END_REF] In this case the region-dependency of the descriptors must be taken into account. If not, 7, 8, 12-15 some correctional terms in the expression of the velocity may be omitted.

Accounting For Region Dependency

In this paper, the region-dependency of descriptors is taken into account in the derivation of the velocity. It is shown that it induces additional terms leading to a greater accuracy of segmentation. The development is general enough to provide a framework for region-based active contour. It is inspired by shape optimization techniques.

J(Ω in ) = α Ω in k in (Ω in , x) dx + β Ωout k out (Ω out , x) dx + Γ k b (x) dx (7)
where α and β are positive constants, Ω in , Ω out , and Γ are such that

   Ω in ∪ Γ ∪ Ω out = image domain D Γ = ∂Ω in = ∂Ω out Ω in = domain inside Γ , (8)
and k in , k out , and k b are positive functions such that

      
k in is minimum in the object, e.g., function (5) k out is minimum in the background, e.g., k out = 0, meaning that no specific information is available for the background k b is minimum on the object boundary, e.g., function

.

Therefore, the minimum of criterion ( 7) is reached if Ω in segments the object of interest. The proposed segmentation method is based on the active contour technique (see Subsection 2.B below). Let us give the notations here. Domains Ω in and Ω out and contour Γ are respectively replaced with dynamic versions depending on evolution parameter τ . At τ equal to zero, an initial domain Ω in (τ = 0) = Ω 0 in is defined, either manually or automatically (equivalently, Ω out (τ = 0) = Ω 0 out and Γ(τ = 0) = Γ 0 ). The active contour converges toward the object boundary as τ increases. It increases by dτ at each iteration when discretizing the evolution equation for computer coding.

B. Differentiating the Criterion

The set of domains Ω in is not a vectorial space. Direct computation of the derivative of criterion (7) with respect to Ω in is not possible. The proposed solution is to use a dynamic scheme where Ω in becomes continuously dependent on an evolution parameter denoted by τ . As τ increases, Ω in (τ ) must act as a minimizing sequence of criterion (7). It is equivalent to finding a family of domain transforms T τ such that T 0 is the identity and

T τ (Ω 0 in ) = Ω in (τ )
where Ω 0 in is an initial contour. As a consequence, the active contour philosophy of the proposed method has a theoretical explanation rather than being used as an implementation tool for an energy minimization problem. Criterion (7) becomes (10) which can be rewritten

J(Ω in (τ )) = α Ω in (τ ) k in (Ω in (τ ), x) dx + β Ωout(τ ) k out (Ω out (τ ), x) dx + Γ(τ ) k b (x) dx
J(τ ) = α Ω in (τ ) k in (τ, x) dx + β Ωout(τ ) k out (τ, x) dx + Γ(τ ) k b (x) dx . (11)
Criterion ( 11) is composed of two types of integrals:

J 1 (τ ) = Ω(τ ) k(τ, x) dx (region integral) (12)
J 2 (τ ) = Γ(τ ) k b (x) dx (contour integral) (13)
For simplicity, k(τ, x) will be written k and k b (x) will be written k b .

Region Integral

Theorem 1 [START_REF] Sokolowski | Introduction to shape optimization: Shape sensitivity analysis[END_REF][START_REF] Delfour | Shapes and geometries: Analysis, Differential Calculus and Optimization[END_REF] Let D be the image domain (Ω(τ ) ⊂ D). Let k be a smooth function on R + × D. Then

dJ 1 dτ (τ ) = J 1 (τ ) = Ω(τ ) ∂k ∂τ dx - ∂Ω(τ ) v • N k dx ( 14
)
where v (actually v(τ, x)) is the velocity of ∂Ω(τ ) and N (actually N (τ, x)) is the inward unit normal to ∂Ω(τ ).

J 1 (τ ) is called the Eulerian derivative of J 1 (Ω(τ )) in the direction of v at τ . It represents the variation of J 1 (τ ) due to both the deformation of integration domain Ω(τ ) according to v and the variation of k. The variation of k is also due to the deformation of domain Ω(τ ).

Corollary 1 [START_REF] Jehan-Besson | DREAM 2 S: Deformable regions driven by an eulerian accurate minimization method for image and video segmentation[END_REF] Let D be the image domain. Let Ω in (τ ), Ω out (τ ), and Γ(τ ) be two domains and a boundary as defined by Eqs. (8). Let k in , respectively k out , be a smooth function on

R + × Ω in (τ ), respectively R + × Ω out (τ ). Let J1 (τ ) be J1 (τ ) = Ω in (τ ) k in (τ, x) dx + Ωout(τ ) k out (τ, x) dx . ( 15
)
Then

d J1 dτ (τ ) = J 1 (τ ) = D ∂K ∂τ dx - Γ(τ ) v • N [[K]] dx - ∂Ωout(τ )\Γ(τ ) w • N Ωout K dx ( 16
)
where

K is the function equal to k in in Ω in (τ ) and k out in Ω out (τ ), [[K]] is the jump of K across Γ(τ ) and it is equal to k in -k out , ∂Ω out (τ )\Γ(τ ) is the boundary of Ω out (τ )
excluding Γ(τ ) (namely ∂D), w is the velocity of ∂Ω out (τ )\Γ(τ ), and N Ωout is the inward unit normal to ∂Ω out (τ )\Γ(τ ).

Corollary 1 is obtained by applying theorem 1 with Ω in (τ ), k in and D = Ω in (τ ) for the first integral of (15) and with Ω out (τ ), k out and D = Ω out (τ ) for the second integral (see Appendix). Note that the image domain being fixed, velocity w is equal to zero and so is the last integral of (16).

Contour Integral

The derivative of J 2 (τ ) is classical.

3 dJ 2 dτ (τ ) = J 2 (τ ) = Γ(τ ) (∇k b • N -k b κ)v • N dx (17)
where κ (actually κ(τ, x)) is the curvature of Γ(τ ).

Complete Criterion

Finally, the derivative of criterion (11) is

J (τ ) = α Ω in (τ ) ∂k in ∂τ dx + β Ωout(τ ) ∂k out ∂τ dx + Γ(τ ) (βk out -αk in + ∇k b • N -k b κ)v • N dx . ( 18
)
The first two integrals in (18) can be reduced to a boundary integral: [START_REF] Jehan-Besson | DREAM 2 S: Deformable regions driven by an eulerian accurate minimization method for image and video segmentation[END_REF] α

Ω in (τ ) ∂k in ∂τ dx + β Ωout(τ ) ∂k out ∂τ dx = Γ(τ ) H(k in , k out )v • N dx (19)
where H(k in , k out ) represents the additional terms mentioned at the beginning of section 2. As a consequence, derivative (18) can be rewritten

J (τ ) = Γ(τ ) ρv • N dx (20) with ρ = H(k in , k out ) + βk out -αk in + ∇k b • N -k b κ . ( 21
)
Velocity v is unknown. It must be chosen so that the value of ( 20) is negative in order to make criterion (11) decrease. Note that this result can be obtained by a classical calculus of variation approach although the development is tedious. [START_REF] Aubert | Image segmentation using active contours: Calculus of variations or shape gradients?[END_REF]

C. Evolution Equation

A way to ensure negativity of derivative ( 20) is to chose v = -ρN . Derivative ( 20) is then

J (τ ) = - Γ(τ ) ρ 2 dx . ( 22
)
Since velocity v can also be written ∂Γ ∂τ , the evolution equation of the active contour is

Γ(0) = Γ 0 ∂Γ ∂τ (τ, x) = -ρ(τ, x)N (τ, x) for all τ ≥ 0 and x ∈ Γ(τ ) . (23)
Equation ( 23) means that, starting from initial contour Γ 0 , active contour Γ evolves until velocity amplitude ρ is equal to zero everywhere along the contour. In other words, active contour Γ at τ + dτ is obtained by deforming active contour Γ at τ according to local velocity -ρN , if it is not equal to zero. If region descriptors k in and k out are region-independent (i.e., if they do not depend on Ω in (τ )), then H(k in , k out ) is equal to zero. However this is not a necessary condition. Note that Eq. ( 23) solves the segmentation problem (the primary goal) and computes the parameters involved in the descriptors simultaneously. These parameters could be used for indexing and retrieval tasks. For example, if k in is defined to be minimum in regions with an homogeneous blue color, the corresponding parameter is the mean value of the blue component within the region. If the segmentation method is applied to a set of RGB images, the parameter values can be stored along with the blue region contours after which the image database can be queried for images containing homogeneous regions with a given blue component value.

D. Examples Of Descriptors

Constant Intensity

An example of a region descriptor that is region-independent is

k in (τ, x) = φ(δ -I(x)) ( 24
)
where φ is a smooth, positive, even function increasing on R + . This descriptor can be used to segment a region of known intensity δ. If k out is defined similarly, then H(k in , k out ) is equal to zero.

Mean Intensity

The mean intensity within Ω in (τ ) is

µ(τ ) = Ω in (τ ) I(x) dx Ω in (τ ) dx . ( 25
)
If descriptor k in is φ(µ(τ ) -I(x)) and k out is equal to zero, then

H(k in , k out )(τ, x) = - µ(τ ) -I(x) Ω in (τ ) dy Ω in (τ ) φ (µ(τ ) -I(y)) dy ( 26
)
If φ is the square function, then H(k in , k out ) is equal to zero. [START_REF] Chan | Active contours without edges[END_REF][START_REF] Debreuve | Space time segmentation using level set active contours applied to myocardial gated SPECT[END_REF] This descriptor can be used to segment an homogeneous region of unknown intensity.

Variance

The variance within Ω in (τ ) is

σ 2 (τ ) = Ω in (τ ) (µ(τ ) -I(x)) 2 dx Ω in (τ ) dx (27)
where µ(τ ) is defined as in (25). If descriptor k in is φ(σ 2 (τ )) and k out is equal to zero, then

H(k in , k out )(τ, x) = φ (σ 2 (τ ))[σ 2 (τ ) -(µ(τ ) -I(x)) 2 ] ( 28
)
Note that H(k in , k out ) depends on x whereas the descriptor (hence, the criterion) does not. This descriptor can be used to segment a region with low noise or low intensity variation.

Shape Of Reference

When segmenting a region using appropriate descriptors, it can be useful to constrain the active contour to stay close to a shape of reference. This can be done using the following criterion 20

J 2 (τ ) = Γ(τ ) k b (τ, x) dx = Γ(τ ) φ(d(x, Γ ref )) dx (29)
where

d(x, Γ ref ) is the signed distance to Γ ref at x, i.e., d(x, Γ ref ) = +|x -y(x)| if x is outside Γ ref -|x -y(x)| otherwise ( 30
)
where y(x) is the closest point to x belonging to Γ ref , and φ is a smooth, positive, even function increasing on R + . In criterion (13), boundary descriptor k b does not depend on τ , i.e., on the boundary. Computation of the derivative of boundary-dependent descriptor ( 29) is different from the developments mentioned in Subsection 2.B. It can be shown [START_REF] Gastaud | Tracking video objects using active contours[END_REF] that the derivative is

J 2 (τ ) = - Γ(τ ) (φ (d)N ref • N + φ(d)κ)v • N dx (31)
where d is a short notation for d(x, Γ ref ) and N ref is the inward unit normal to Γ ref at y(x). The deformation between the shape of reference and the segmentation contour is a free-form deformation (i.e., deformation with unlimited degrees of freedom) as opposed to a parametric transform [START_REF] Chen | Using prior shapes in geometric active contours in a variational framework[END_REF][START_REF] Cremers | Diffusion snakes: Introducing statistical shape knowledge into the mumford-shah functional[END_REF] (i.e., global deformation with few degrees of freedom, e.g., a combination of a translation, a rotation, and a scaling). Although the derivatives of the criteria may be similar in both cases, the parametric transform approach tends to turn the segmentation problem into an iterative best fit search among the set of contours obtained by parametric transform.

Active Contour Implementation

A. From Parametric To Implicit

The first implementations of the active contour technique were based on a parametric, or explicit, description of the contour. [START_REF] Kass | Snakes: Active contour models[END_REF] It corresponds to a Lagrangian approach. However, management of the evolution, particularly topology changes and sampling density along the contour, is not simple. [START_REF] Charbonnier | Une étude des contours actifs : modèles classique, géométrique et géodésique[END_REF] Instead, an Eulerian approach known as the level set technique [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations[END_REF][START_REF] Barles | Remarks on a flame propagation model[END_REF] can be used. In two dimensions, the contour is implicitly represented as the intersection of a surface u = f (x), where x ∈ D and elevation f is a continuous function with positive and negative values, with the plane of elevation zero.

The contour can also be seen as the isocontour of level zero on the surface. In three dimensions, the contour is the isosurface of level zero in a volume. In n dimensions, the contour is the hyperplane of level zero with the space filled in with the values of a real, continuous function. Note that the contour can actually be composed of several closed contours without intersections with each other. By a continuous change of function f , a contour can appear or disappear without explicit handling. Unfortunately, even with the narrow band implementation, [START_REF] Chopp | Computing minimal surfaces via level set curvature flow[END_REF][START_REF] Adalsteinsson | A fast level set method for propagating interfaces[END_REF] the level set technique has a rather high computational cost and extension of the velocity to levels other than the zero level is not straightforward [START_REF] Adalsteinsson | The fast construction of extension velocities in level set methods[END_REF][START_REF] Gomes | Reconciling distance functions and level sets[END_REF] while it is theoretically necessary. Moreover, a curvature term (minimum length penalty) is usually added to the velocity expression in order to decrease the influence of image noise on the evolution. However, the curvature being a second derivative term, its numerical approximation is usually not accurate. Similarly, the level set gradient (equivalent to the isocontour normal) is numerically approximated. The fast marching technique [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF] is another Eulerian approach to the active contour representation. It could be seen as the limit case where the narrow band is one grid element thick. This technique is much faster than the level set technique. However it can be applied only if one can guarantee that the velocity expression keeps the same sign during evolution.

B. Splines: Back To The Parametric Approach

A cubic B-spline has several interesting properties: It is a C 2 curve, [START_REF] Thevenaz | Interpolation revisited [medical images application[END_REF] it is an interpolation curve that minimizes a term close to the square curvature, [START_REF] Jacob | A unifying approach and interface for splinebased snakes[END_REF] and it has an analytical equation (depending on control points) between each pair of consecutive sampling points (p i , p i+1 ). It minimizes the following criterion with the constraint that it passes through the sampling points

J CubicS (Γ) = Γ Γ"(x) 2 dx ( 32
)
where Γ" is the second derivative of Γ with respect to the arc length parameter. Note that the sampling points correspond to the ends, or knots, of the spline segments. The velocity has to be computed at the sampling points only. If the sampling is regular, the normal and the curvature can be computed using an exact, fast, recursive filtering algorithm applied to the control points. Therefore, the spline implementation is much less time consuming than the level set technique. [START_REF] Precioso | B-spline active contours with handling of topological changes for fast video segmentation[END_REF][START_REF] Precioso | Regular spatial B-spline active contours for fast video segmentation[END_REF] Moreover, the minimum curvature-type term property helps in decreasing the influence of noise without the need to add a curvature term to the velocity. Nevertheless, noise in the image still implies noise in the velocity which, if sampling is fine, usually leads to an irregular contour because, despite the smooth curvature property, a cubic B-spline is an interpolation curve. Finally, let us recall that, unlike with level sets, management of topology changes has to be performed explicitly and is not an easy task. Also note that the number of sampling points is a parameter. Less sampling points implies a curve more rigid, or smoother, preventing from segmenting regions with a high curvature.

C. Smoothing Splines

A smoothing spline is an approximation curve controlled by a parameter balancing the trade-off between interpolation error and smoothness. [START_REF] Unser | B-spline signal processing: Part I-theory[END_REF][START_REF] Precioso | Smoothing B-spline active contour for fast and robust image and video segmentation[END_REF] It minimizes the following criterion

J SmoothS (Γ) = λ Γ Γ"(x) 2 dx + i (p i -Γ i ) 2 (33)
where λ is a parameter, p i is a sampling point, and Γ i is a knot. Note that knots do not correspond to sampling points anymore. The smoothing spline is only constrained to be confined to a "band" surrounding the sampling points. The smaller parameter λ, the narrower the band. When λ is equal to zero, the smoothing spline is a classical interpolation spline. Thus, if sampling is fine and noise is high, the smoothing spline can still be smooth. As with cubic B-splines, normal and curvature can be computed exactly and efficiently.

Examples Of Applications

The following results were obtained using the level set technique except otherwise noted.

A. Image Segmentation

Region Competition

Image "X-ray" is a 512 × 512, greyscale image. Descriptor k in and k out were both variance descriptors (see Subsection 2.D.3). This type of segmentation is called a region competition. To give a brief comparison in terms of computation time between the level set technique and the smoothing spline implementation, the image was segmented using both algorithms. See Fig. 1. The segmentation using level sets was Fig. 1. Segmentation of "X-ray": Segmentation using a smoothing spline with 128 sampling points (top left), with 256 sampling points (top middle), and with 512 sampling points (top right) and λ equal to 0.005; Segmentation using level sets (bottom).

performed roughly 15 times slower than the segmentation using a smoothing spline with 256 sampling points. With 128 sampling points the ratio increases to 22 times. However the segmentation contour is less accurate because it is too smooth compared to the shape of the object of interest. With 512 sampling points the ratio decreases to 6.5 times and the segmentation contour is not smooth enough. The influence of parameter λ can be seen on Fig. 2. The segmentation contour is smoother with a higher λ.

Fig. 2. Segmentation of "X-ray": Segmentation using a smoothing spline with 256 sampling points and λ equal to: 0.005 (left) and 0.025 (right).

Shape Of Reference

Image "Phone" is a 352 × 288, color image. Two descriptors were used for the segmentation: A variance descriptor (see Subsection 2.D.3) for the face (k in /statistical descriptor) and a shape of reference descriptor (see Subsection 2.D.4) as a constraint (k b /geometrical descriptor). Descriptor k out was taken equal to zero. Combination of these descriptors implies a competition between the shape prior and the statistical information of the object to be segmented. If the shape of reference constraint is omitted, the face segmentation includes part of the hand of the character and does not include the lips. A shape of reference was heuristically defined allowing to segment accurately the face. See Fig. 3.

B. Sequence Segmentation

Sequence "Akiyo" is composed of three hundred 352 × 288, color images. Descriptor k in was taken equal to a constant penalty. Descriptor k out was defined as the difference between the current image and a robust estimate of the background image computed using the n b previous images [START_REF] Jehan-Besson | Video object segmentation using eulerian region-based active contours[END_REF] (n b = 20), or less for the first n b images. This descriptor takes advantage of the temporal information of the sequence: It is a motion detector.

In case of a moving background due to camera motion, a mosaicking technique can be used to estimate the background image. [START_REF] Gastaud | Video segmentation using region based active contours on a group of pictures[END_REF] Descriptor k out was taken equal to zero. The first image of the sequence was not segmented since no previous images were available to estimate the background. See Fig. 4 and5.

C. Tracking

Sequence "Erik" is composed of fifty 256 × 256, color images. Given a segmentation of Erik's face in the first image, the purpose was to track the face throughout the sequence using the segmentation of the previous image to constrain the segmentation of the current frame. [START_REF] Gastaud | Tracking video objects using active contours[END_REF] See Fig. 6. Two descriptors were used: A variance descriptor (see Subsection 2.D.3) for the face (k in /statistical descriptor) and a shape of reference descriptor (see Subsection 2.D.4) as a constraint (k b /geometrical descriptor). Descriptor k out was taken equal to zero. The shape of reference in the current image was defined as an affine transform (i.e., a combination of a translation, a rotation, and a scaling) of the segmentation contour in the previous image. This transform can be interpreted as the global motion of the object of interest between the previous and

Fig. 3 .

 3 Fig. 3. Segmentation of "Phone": Segmentation without shape of reference (left); Shape of reference (middle); Segmentation with the shape of reference constraint (right).

Fig. 4 .

 4 Fig. 4. Segmentation of Akiyo: Image 2, evolution from the initial contour until convergence.

Fig. 5 .

 5 Fig. 5. Segmentation of Akiyo: Images 5, 15, 25, 40, 55, and 70 out of 300.

Fig. 6 .vvv

 6 Fig. 6. Tracking of Erik's face: Images 1, 10, 27, 40, and 50 out of 50.

  [START_REF] Jehan-Besson | DREAM 2 S: Deformable regions driven by an eulerian accurate minimization method for image and video segmentation[END_REF][START_REF] Gastaud | Tracking video objects using active contours[END_REF]

Acknowledgments

Image "Phone" was provided by france telecom, France. The X-ray image was provided by Diagnostic Medical Systems (DMS), Montpellier, France.

Appendix

Applying theorem 1 with Ω in (τ ), k in and D = Ω in (τ ) for the first integral of (15) and with Ω out (τ ), k out and D = Ω out (τ ) for the second integral, we have

Boundary ∂Ω out (τ ) is the union of ∂D and Γ(τ ). Note that the velocity of Γ(τ ) is identical whether it is seen from Ω in (τ ) or Ω out (τ ). On the other hand, the normal is