M L Delle Monache
email: ml.dellemonache@inria.fr

T Liard

A Rat

R Stern

R Badhani

B Seibold
email: sei-bold@temple.edu

J Sprinkle
email: sprinkle@ece.arizona.edu

D Work
email: dan.work@vanderbilt.edu

B Piccoli
email: piccoli@camden.rutgers.edu

Feedback control algorithms for the dissipation of traffic waves with autonomous vehicles

This article considers the problem of traffic control in which an autonomous vehicle is used to regulate human piloted traffic to dissipate stop and go traffic waves. We first investigate the controllability of well-known microscopic traffic flow models, namely i) the Bando model (also known as the optimal velocity model), ii) the follow-the-leader model, and iii) a combined optimal velocity follow the leader model. Based on the controllability results, we propose three control strategies for an autonomous vehicle to stabilize the other, human-piloted traffic. We subsequently simulate the control effects on the microscopic models of human drivers in numerical experiments to quantify the potential benefits of the controllers. Based on the simulations, finally we conduct a field experiment with 22 human drivers and a fully autonomous-capable vehicle, to assess the feasibility of autonomous vehicle based traffic control on real human piloted traffic. We show that both in simulation and in the field test that an autonomous vehicle is able to dampen waves generated by 22 cars, and that as a consequence, the total fuel consumption of all vehicles is reduced by up to 20%.

Introduction

Currently, the vehicular transportation system is undergoing a major transition from vehicles in which humans are responsible for all driving tasks, to one in which automation is responsible for all driving tasks. The transition is defined [START_REF]Taxonomy and definitions for terms related to on-road motor vehicle automated driving systems[END_REF] in terms of various levels of automation. The levels range from level one autonomous vehicles (AVs) available today that provide the driver with minor technological assistance (e.g., stability control or lane correction assist), to level five AVs which operate autonomously in all scenarios and in which humans cannot intervene. As the penetration rate of vehicles at each level of automation increases and shifts up the scale, new opportunities are arising to use automated vehicles to begin controlling the overall traffic flow.

A paradigm of traffic control in which some automated vehicles are also acting as traffic control devices is beginning to emerge. In particular, the works [START_REF] Davis | Effect of adaptive cruise control systems on traffic flow[END_REF][START_REF] Talebpour | Influence of connected and autonomous vehicles on traffic flow stability and throughput[END_REF][START_REF] Guériau | How to assess the benefits of connected vehicles? a simulation framework for the design of cooperative traffic management strategies[END_REF][START_REF] Wang | Cooperative car-following control: Distributed algorithm and impact on moving jam features[END_REF][START_REF] Knorr | Reducing traffic jams via VANETs[END_REF], explore the possibility of adaptive cruise controlled vehicles (e.g., level one automation) to influence traffic flow, for example by smoothing the flow and/or increasing the flow rate. On the extreme end, when the adaptive cruise controlled vehicles are also endowed with communication capabilities between vehicles (referred to as Cooperative Adaptive Cruise Control (CACC) systems), small headways can be achieved and substantial increases in freeway throughput can be obtained [START_REF] Levine | On the optimal error regulation of a string of moving vehicles[END_REF][START_REF] Swaroop | String stability of interconnected systems[END_REF][START_REF] Shladover | Review of the state of development of advanced vehicle control systems (avcs)[END_REF][START_REF] Fenton | Automated highway studies at the Ohio State University-an overview[END_REF][START_REF] Darbha | Intelligent cruise control systems and traffic flow stability[END_REF][START_REF] Besselink | String stability and a delay-based spacing policy for vehicle platoons subject to disturbances[END_REF][START_REF] Ioannou | Intelligent cruise control: theory and experiment[END_REF][START_REF] Buehler | The DARPA urban challenge: autonomous vehicles in city traffic[END_REF][START_REF] Rajamani | Design and experimental implementation of control for a platoon of automated vehicles[END_REF]. On the experimental side, field experiments with commercial adaptive cruise control vehicles illustrate current technology may in fact amplify traffic oscillations [START_REF] Milanés | Cooperative adaptive cruise control in real traffic situations[END_REF], while correctly designed CACC systems dissipate these oscillations [START_REF] Milanés | Cooperative adaptive cruise control in real traffic situations[END_REF][START_REF] Ge | Experimental validation of connected automated vehicle design among human-driven vehicles[END_REF]. An experiment to harmonize speeds on a US freeway was recently reported [START_REF] Ma | Freeway speed harmonization[END_REF][START_REF] Lu | Novel algorithm for variable speed limits and advisories for a freeway corridor with multiple bottlenecks[END_REF].

The use of vehicles as traffic controllers may also integrate with more classical traffic control infrastructure, such as ramp meters and variable speed limit systems [START_REF] Nissan | Evaluation of the impact of advisory variable speed limits on motorway capacity and level of service[END_REF][START_REF] Hegyi | Optimal coordination of variable speed limits to suppress shock waves[END_REF][START_REF] Smulders | Control of freeway traffic flow by variable speed signs[END_REF][START_REF] Hegyi | SPECIALIST: A dynamic speed limit control algorithm based on shock wave theory[END_REF][START_REF] Popov | Distributed controller design approach to dynamic speed limit control against shockwaves on freeways[END_REF][START_REF] Papageorgiou | Freeway ramp metering: An overview[END_REF][START_REF] Gomes | Optimal freeway ramp metering using the asymmetric cell transmission model[END_REF][START_REF] Papageorgiou | ALINEA: A local feedback control law for on-ramp metering[END_REF], or systems which combine the two strategies [START_REF] Hegyi | Model predictive control for optimal coordination of ramp metering and variable speed limits[END_REF][START_REF] Papamichail | Integrated ramp metering and variable speed limit control of motorway traffic flow[END_REF][START_REF] Lu | Combining variable speed limits with ramp metering for freeway traffic control[END_REF][START_REF] Han | Resolving freeway jam waves by discrete first-order model-based predictive control of variable speed limits[END_REF]. One limitation of the infrastructure based solutions is their limited spatial resolution, as well as the need for driver compliance in the case of speed advisory based systems. The interest to use automated vehicles for traffic control is partly motivated by these limitations.

Considering the new direction of traffic control in which automated vehicles act as actuators, the main contribution of this article is to propose and assess control algorithms designed to dissipate stop and go traffic waves with an autonomous vehicle. Stop and go traffic waves are present on freeways and have many triggering events, such as lane changes. Strikingly, in the seminal experiments of Sugiyama et al. [START_REF] Sugiyama | Traffic jams without bottlenecks -experimental evidence for the physical mechanism of the formation of a jam[END_REF] and Tadaki et al. [START_REF] Tadaki | Phase transition in traffic jam experiment on a circuit[END_REF], human driving behavior alone was shown to be sufficient to trigger stop and go waves. The experiments were conducted on a single lane circular track with real human drivers, and the uniform flow at the experiment start quickly breaks down into a persistent stop and go wave that travels against the flow of traffic.

Motivated by this experiment, we design control algorithms to be implemented on an autonomous vehicle with the goal of dissipating stop and go waves caused by human driving behavior. We proceed as follows. First we model the vehicular traffic at the microscopic scale using one or a combination of two well-known microscopic models, namely i) the Bando or optimal velocity (OV) model [START_REF] Bando | Dynamical model of traffic congestion and numerical simulation[END_REF] and the ii) follow the leader (FTL) [START_REF] Reuschel | Vehicle movements in a platoon with uniform acceleration or deceleration of the lead vehicle[END_REF][START_REF] Reuschel | Vehicle movements in a platoon[END_REF][START_REF] Pipes | An operational analysis of traffic dynamics[END_REF]. With the models defined, we show that traffic described by a linerization of the FTL model is not controllable via a single AV, implying the nonlinear FTL model is not linearly controllable. This result prevents the use of simple linear controllers to stabilize the traffic around the uniform speed equilibrium traffic state. On the other hand, we show that traffic described by the optimal velocity model is locally controllable by an autonomous vehicle. We show that for driving dynamics that include both the optimal velocity and follow the leader terms, the resulting model is also locally controllable 1We then proceed to design three different controls in which the AV is used to dampen stop and go waves. The first two are based on Lyapunov functions and only require measurements of the AV speed and the speed of nearby vehicles, while the third is a PID-type control in which the AV is controlled using only measurements of its own speed over several proceeding timesteps.

We first assess the effectiveness of the wave dampening controllers though numerical experiments. In addition to showing a reduction in the wave strength, we also quantify the benefits in terms of a reduction of total fuel consumed by all vehicles in the simulation. Based on the positive performance in simulation, we proceed to field validate one control algorithm with a dive-by-wire autonomous capable vehicle on a track with 22 vehicles driven by humans. The experiment shows that stop and go waves can be dampened, and the projected reduction of fuel consumption from the simulations (approximately 20%) is confirmed via real-time fuel consumption loggers installed on the experimental vehicles.

The remainder of the article is organized as follows. In Section 2, we review the main microscopic traffic flow models investigated in this work. Section 3 establishes the main controllability results of the models, while Section 4 describes the design and testing the wave dissipating controllers in simulation. In Section 5 we test one of the wave dampening controllers in a field experiment with real human drivers and an autonomous vehicle. Limitations and future directions are explored in Section 6.

Microscopic traffic models

Traffic models are usually defined in categories, depending on the scale at which they represent vehicular traffic, including microscopic, mesoscopic, macroscopic and cellular. For a review of models at various scales see [START_REF] Helbing | Traffic and related self-driven many-particle systems[END_REF][START_REF] Bellomo | On the modeling of traffic and crowds: A survey of models, speculations, and perspectives[END_REF][START_REF] Garavello | Models for vehicular traffic on networks[END_REF][START_REF] Piccoli | Vehicular Traffic: A Review of Continuum Mathematical Models[END_REF]. Moreover, some approaches are based on model-agnostic simulation tools, such as deep learning and neural networks [START_REF] Hongfei | Develop a car-following model using data collected by" five-wheel system[END_REF][START_REF] Chong | Simulation of driver behavior with agent-based backpropagation neural network[END_REF][START_REF] Khodayari | A modified car-following model based on a neural network model of the human driver effects[END_REF][START_REF] Zhou | A recurrent neural network based microscopic car following model to predict traffic oscillation[END_REF].

In this paper we focus on microscopic models. Microscopic models are suitable for in-silico verification before experimental testing since they describe human driving behavior at the individual vehicle level. One such model is the combined Bando and follow the leader model, which can be formulated as ẋi = vi, vi = f (xi+1, xi, vi+1, vi), where xi is the position of i-th car, vi its velocity, and i + 1 is the index of the car ahead. More precisely, for what concerns the follow the leader model, it was introduced in [START_REF] Reuschel | Vehicle movements in a platoon with uniform acceleration or deceleration of the lead vehicle[END_REF][START_REF] Reuschel | Vehicle movements in a platoon[END_REF][START_REF] Pipes | An operational analysis of traffic dynamics[END_REF] and it assumes that the acceleration of a vehicle is given by the neighbouring vehicles. The main influence comes from the next vehicle, whose index is i + 1, that is also called leading vehicle. The main dynamics described by this model is given by:

   ẋi = vi vi = vi+1 -vi xi+1 -xi . 1 i N (1)
This model has the following properties:

• The acceleration depends on the relative velocity ∆v = vi+1 -vi

• The velocity vi(t) of the vehicle depends on the velocity of the vehicle in front such that the distance from the vehicle in front is safe.

A drawback of these models is that the acceleration is zero when the relative velocity is zero independently of the headway d ≡ xi+1 -xi That is, extremely small headways are allowed even when travelling with extremely high speed, [START_REF] Orosz | Traffic jams: dynamics and control[END_REF]. A model that fixed this problem is the optimal velocity model (introduced by [START_REF] Bando | Dynamical model of traffic congestion and numerical simulation[END_REF]) that describes the adaptation of the actual speed to the optimal velocity V (•) which stands for the desired speed defined by

V (x) = vmax tanh(x -lv -ds) + tanh(lv + ds) 1 + tanh(lv + ds) , (2)
where lv is the length of cars and ds > 0 is the safe distance between cars. The optimal velocity has the property that it tends to zero for small headways and it achieves the maximum value for large headways.

The full model reads:

ẋi = vi vi = V (xi+1 -xi) -vi. 1 i N (3)
In this model, a driver controls acceleration or deceleration according to the difference between the optimal velocity and his own velocity. The equilibrium point for this model is achieved when all cars at constant speed and with the same headway, [START_REF] Konishi | Decentralized delayed-feedback control of an optimal velocity traffic model[END_REF]. For this model it is possible to derive a feedback law such that the controlled traffic system is stable, [START_REF] Nakayama | Effect of looking at the car that follows in an optimal velocity model of traffic flow[END_REF].

Controllability results for microscopic models

In this Section, we provide theoretical results about controllability for the microscopic models introduced in Section 2. More precisely, we will focus on a ring-road setting, which reproduce the situation of the celebrated Sugiyama experiment [START_REF] Sugiyama | Traffic jams without bottlenecks-experimental evidence for the physical mechanism of the formation of a jam[END_REF] with a single AV, which can be controlled, and investigate the controllability of the corresponding control system. Let us first recall some basic facts about control systems, referring the reader to the books [START_REF] Bressan | Introduction to the Mathematical Theory of Control[END_REF] and [START_REF] Coron | Control and nonlinearity[END_REF] for details. A control systems is a dynamical system written as:

ẏ = f (y, u) (4)
where y ∈ IR n represents the state of the system and u ∈ U ⊂ IR m represents the control vector, i.e. the parameters which can be chosen by an external agent. In our setting y represent the state of the traffic model (e.g. position and velocity of cars) and u the acceleration or speed of the AV which can be controlled.

A system (4) is said controllable if for every states y1, y2, there exists T > 0 and a control function ū : [0, T ] → U such that the solution to the Chauchy problem ẏ = f (t, ū(t)), y(0) = y1 satisfies y(T ) = y2. In other words we can steer the system from y1 to y2 in time T > 0 with a suitable control ū. Similarly a system is said locally controllable at ȳ if for every δ > 0 sufficiently small there exists T > 0 such that we can steer the system from ȳ to any y with |y -ȳ| < δ in time T . In other words we can reach sufficienly close states in uniformly bounded time.

n 1 2 3 n -1 n -2 v 2 x 2 ≈ l v
An important tool to investigate controllability (and local ones) is linearization. Given a control system (4), such that U contains a neighborhood of 0, and ȳ, we can consider the linearized system at ȳ:

ż = A • z + B • u, (5)
where A = Dyf (ȳ, 0) (the Jacobian matrix of f w.r.t. y computed at (ȳ, 0)) and B = Duf (ȳ, 0) (the Jacobian matrix of f w.r.t. u computed at (ȳ, 0)). For linear systems, there is a simple criterion for controllability. First, given the linear control system (5), define the Kalman controllability matrix :

K(A, B) = [B, A • B, . . . , A n-1 • B]. ( 6
)
Since A is an n × n matrix and B is a n × m matrix, then K(A, B) is an n × n • m matrix. We can now state the follwing (see [4, Theorem 3.6.2]):

Theorem 1. (Kalman controllability theorem) The system (5) is controllable if and only if the matrix K(A, B) has full rank (i.e. equal to n).

The fact that we can limit ourselves to the exponent n -1 in the definition of K(A, B) follows from Cayley-Hamilton Theorem, which ensures that A is root of its characteristic polynomial, i.e. there exists αi such that A n = i A i = 0 n-1 αiA i . The local controllability of a system can be investigated by looking at its linearization. More precisely we have the following (see [START_REF] Bressan | Introduction to the Mathematical Theory of Control[END_REF]Theorem 3.7.1]):

Theorem 2. Consider a control system (4) and its linearization (5) at a given state ȳ. If (5) is controllable then (4) is locally controllable at ȳ. The converse of this Theorem is not true: a nonlinear control system may be controllable even if its linearization fails to be controllable. However, the controllability of the linearization ensures the existence of linear stabilizing feedbacks, i.e. control laws of the type u(y) = K • y that renders the system (4) locally asymptotically stable at ȳ, see [START_REF] Bressan | Introduction to the Mathematical Theory of Control[END_REF]Theorem 4.2.3].

We are now ready to define our control system for traffic control on a ring-road via one AV, see Figure 1. Let (a, b) ∈ IR 2 \{(0, 0)} and consider the control system of n vehicles along a ring-road of length L described by the FTL-Bando model:

         ẋi = vi, 1 i n, vi = a vi+1 -vi (xi+1 -xi) 2 + b [V (xi+1 -xi) -vi] , 1 i n -1, vn = u, (7)
where xi is the position of the i th vehicle, vi its velocity and V (•) stands for the desired speed which is defined in 2. Notice that we can control only the acceleration u of the n-th vehicle, which represents the AV.

Remark 1. If b = 0, the nonlinear control system ( 7) is the Follow-The-Leader model (FTL). If a = 0, the system ( 7) is the optimal velocity model.

Our aim is to steer the system, ( 7) to a speed equilibrium state, i.e. a state so that all vehicles have the same velocity v. If b = 0, i.e. if there is a non vanishing optimal velocity term, a speed equilibrium corresponds to all vehicles which are not controlled (i.e., all vehicles but the n-th ones) having the same headway distance d, while the AV may have a different headway. Moreover, to be in speed equilibrium we must have V (d) = v, thus d is fixed once v is fixed. On the other side if b = 0 then all states with same speed are an equilibrium with any set of headways. We summarize this analysis in the following Lemma:

Lemma 1. Consider the control system (7) and a fixed speed v. If b = 0 then the following holds. Let d be such that

V (d) = v. If d L n then the only speed equilibrium is given by x(t) such that xi+1(t) -xi(t) ≡ d, i = 1, . . . , n -1 and ẋi(t) ≡ v, i = 1, . . . , n. Otherwise there exists no speed equilibrium. If b = 0 then every x(t) such that ẋi(t) ≡ v, i = 1, . . . , n, is an equilibrium. Now fix (d, v) ∈ IR 2 + , with V (d) = v if b = 0, then we rewrite (7) as Ẏ = f (Y, u). (8)
where

Y = (yi)i=1,••• ,2n-1 = (x2 -x1 -d, . . . , xn -xn-1 -d, v2 -v1, . . . , vn -v) and f is defined by f (Y, u) =                    yn . . . y2n-2 a y n+1 (y 2 +d) 2 -yn (y 1 +d) 2 + b [V (y2 + d) -V (y1 + d) -yn] . . . a y 2n-2 (yn-1+d) 2 - y 2n-3 (yn-2+d) 2 + b [V (yn-1 + d) -V (yn-2 + d) -y2n-3] u -a y 2n-2 (yn-1+d) 2 -b [V (yn-1 + d)-( y2n-1 -y2n-2 + v )] u                    (9)
Notice that the vector Y has 2n -1 component: n -1 differences between the headway of the n -1 not controlled vehicles and the equilibrium headway d, n -1 velocity differences with leading vehicle for the not controlled vehicles and the difference of AV velocity with the equilibrium ones. Moreover we have the relation x1 -xn = L -n-1 i=1 yi + (n -1)d expressing the headway of the AV. We are interested in speed equilibria with equispaced cars, thus we notice that Y * = 0 ∈ IR 2n-1 is the only equilibrium of (8) if b = 0 and u = 0. In other words (Y * , 0) is a speed equilibrium of the controlled system (8), i.e. f (Y * , 0) = 0. Notice that for d = L n , V (d) = v, at the speed equilibrium (0, 0) all vehicles are equispaced and drive at the same speed vi = v.

We are now ready to state the controllability results for the linearized systems at equilibria. First we have the following:

Theorem 3. Let b = 0, a = 0 and (d, v) ∈ IR 2
+ . The linearization of system (8) at the speed equilibrium (Y * , 0) ∈ IR 2n is not controllable.

In other words the linearization FTL model at speed equilibrium is not controllable. This does not prevent the nonlinear system to be controllable, but prevents the use of simple linear controls to drive the system to equilibrium.

Our next result is the following:

Theorem 4. Let a = 0, b = 0, (d, v) ∈ IR 2 + with V (d) = v.
The linearization of system (8) at the speed equilibrium (Y * , 0) ∈ IR 2n is controllable. Therefore, the nonlinear control system (8) is locally controllable at (Y * , 0).

Finally for the combined FTL-OV model we have the following:

Theorem 5. We assume that n 9. Let (a, b) = 0, (d, v) ∈ IR 2 + with V (d) = v. If a = (y * + d)V (y * + d
) and b = 0, the linearization of system (8) at the equilibrium point (Y * , 0) ∈ IR 2n is controllable. Therefore, the nonlinear control system (8) is locally controllable at (Y * , 0). Remark 2. Since the Optimal Velocity model is locally controllable, one expect the same to be true for the combined FTL-OV model, except possibly some resononant value of the parameters. This is exactly what it proved in Theorem 5. The complexity of the system does not allow to deal with arbitrary large dimensione (only for n 9), but we conjecture that, for n > 9, the nonlinear control system (8) is still locally controllable at (Y * , 0). The proofs of the Theorems are postponed to the Appendix.

Controls and simulations

In this Section we describe control algorithms to stabilize a FTL-Bando model to a speed equilibrium and test them in silico, i.e. via simulations, on a setting reproducing the Sugiyama experiment [START_REF] Sugiyama | Traffic jams without bottlenecks-experimental evidence for the physical mechanism of the formation of a jam[END_REF]. We start defining a Proportional-Integral-Derivative (PID)-type control, which is based on velocity measures of the AV over a fixed time horizon. The control includes saturation terms to avoid collisions and too large headways for the AV. Then we define controls based on Lyapunov functions, using the fact that the system is control affine, see see [8, Definition 3.12], i.e. the control appear linearly with a vector coefficient depending on state. Such techniques are usually referred to as Jurdjevic-Quinn controls [START_REF] Jurdjevic | Controllability and stability[END_REF].

PID control

The idea behind this controller is that the autonomous vehicle may estimate the average speed of the vehicles in front, and then drive according to the average speed, safety permitting. An estimate of the average speed required by the controller is obtained by measuring the autonomous vehicle speed over a large enough time horizon. Note that this requires that there are several waves present, not just a single one, so that the past is informative of the future.

The controller determines a command velocity u following a standard proportional integral control logic. In order for the controller to be efficient it needs to be augmented with saturation: for small gaps the autonomous vehicle should follow the lead vehicle speed to avoid dangerous situations, while for large gaps, the autonomous should catch up to the lead vehicle.

More precisely, this controller estimates the desired velocity, V d , as a temporal average of the autonomous vehicle's own velocity over an interval. Letting v AV 1 , . . . , v AV m denote the autonomous vehicles velocities over the last m measurements, the desired velocity is computed as the temporal average

V d = 1 m m j=1 v AV j .
In practice, we choose m corresponding to a 38 second interval, which is approximately the time required to travel one lap around the ring.

The desired average velocity is then translated into a target velocity depending on the current gap between the autonomous vehicle and lead vehicle:

v target = V d + 1 m s × min(max( ∆x-7m 23m , 0), 1) , (10)
This allows the autonomous vehicle to drive faster than the average velocity and catch up to the lead vehicle, should it face a big gap, while at lower gaps the target velocity reduces to the average V d . The commanded velocity is updated via:

uj+1 = βj(αjv target j + (1 -αj)v lead j ) + (1 -βj)uj , (11)
where the subscript j denotes the time step. This rule [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF] chooses the new commanded velocity as a weighted average of the prior commanded velocity, the target velocity, and the lead vehicle's velocity.

The weights αj and βj depend on the gap as follows:

αj = min(max( ∆x-∆x s γ , 0), 1) (12)
In [START_REF] Fenton | Automated highway studies at the Ohio State University-an overview[END_REF], the distance ∆x s is a safety distance. We have αj = 0 if ∆x ≤ ∆x s and αj = 1 if ∆x ≥ ∆x s + γ, meaning that for relatively short gaps, only the lead vehicle's velocity matters, while for relatively large gaps, only the target velocity is averaged with the commanded velocity. The parameter γ controls the the rate at which α transitions from 0 to 1, and is set to γ = 2 m in the current implementation. This means that when the gap is short, the autonomous vehicle has the same speed of the lead vehicle, while when the gap is larger the autonomous speed tends towards the target vehicle, which allows the autonomous vehicle to reduce the gap with the lead vehicle. The parameter βj determines how rapidly the controller adjusts to new situations (with more rapid adjustments occurring in more safety-critical situations). At its core, this is a PID controller, but with a saturation at small gaps (for safety purposes), and a saturation at large gaps (so that the autonomous vehicle closes gaps).

The safety distance is implemented as ∆x s = max(2 s × ∆v, 4 m). The term 2 s × ∆v represents the recommended safe following headway of 2 s, with a lower bound of 4 m.

Lyapunov functions and controls

The system ( 8) is a control affine system (see [START_REF] Coron | Control and nonlinearity[END_REF]Definition 3.12]). More precisely, (8) can be written as

Ẏ = f a,b (Y ) + uf1(Y ),
with f a,b (Y ) given by ( 9) with u = 0 and and

f1(Y ) =        0 . . . 0 1 1        .
Since we want to steers the system towards the speed equilibrium, we can consider the control Lyapunov function V1(Y ) = vn -v = y2n-1. Then its derivative along trajectories is given by ∇V1

• (f a,b (Y ) + uf1(Y )) = e2n-1 • (f a,b (Y ) + uf1(Y )) = u.
To maximize the decrease we can simply choose

u1 = -y2n-1 = -(vn -v). ( 13
)
Such control is smooth, vanishes at 0, and globally asymptotically stabilizes to the speed equilibrium.

Remark 3. To have a different gain we may choose the control uα = -αy2n-1 = -α(vn -v) with α > 0.

Let us now focus on last two components:

y2n-2 y2n-1 = u -a y 2n-2 (y * +d) 2 u . (14)
We introduce the control Lyapunov function:

V2(Y ) = y 2 2n-1 2 + y 2 2n-2 2 , for every Y = y2n-2 y2n-1 ∈ IR 2 .
For every Y ∈ IR 2 \{(0, 0)}, choosing u = -sign(y2n-1 + y2n-2) with > 0, we have,

∇V2(Y ) • f (Y, u) = (y2n-2 + y2n-1)u -a y 2 2n-2
(y * +d) 2 < 0. Thus, V2 is a control Lyapunov function for the control system (14) satisfying the small control property (see [8, definition 12.1]). Using [START_REF] Coron | Control and nonlinearity[END_REF]Theorem 12.4]), the control

u2 = -(y2n-1 + y2n-2) = -vn - vn-1 + v 2 , (15)
is smooth, vanishes at 0 ∈ IR 2 and globally asymptotically stabilizes the control system [START_REF] Ge | Experimental validation of connected automated vehicle design among human-driven vehicles[END_REF].

Simulations

In this section we demonstrate the capabilities of the control laws described in the previous sections via numerical simulations. The parameters used for the simulations are as follows. We consider N = 22 vehicles, with one autonomous vehicle and N = 21 human-driven vehicles. We consider that the humandriven vehicles follow the dynamics described in [START_REF] Chong | Simulation of driver behavior with agent-based backpropagation neural network[END_REF] with the following parameters a = 0.5m/s 2 and b = 20s -1 . The speed vmax = 9.75m/s and the vehicle length is chosen to be lv = 4.5m. Such parameters allow to fit the data from the experiment of [START_REF] Sugiyama | Traffic jams without bottlenecks-experimental evidence for the physical mechanism of the formation of a jam[END_REF], see Figure 2. In the Figure, the red trajectory corresponds to the autonomous vehicle and the grey ones to human-driven vehicles. One can notice the appearance of strong stop and go waves that start at time 60s and are propagated on the ring. Visually, there have the effect of "wrinkles" in the pattern.

For comparison, we can now see the effects of the different controls. For the autonomous vehicle we will chose a different dynamics according to the control strategy that we are going to simulate. The simulation begins with all vehicles using the human driving dynamics described in [START_REF] Chong | Simulation of driver behavior with agent-based backpropagation neural network[END_REF] with no control. Control on the AV is activated after t = 40s.

In Figure 3, we can see the effects of the PID control on vehicle trajectories and on the speed profile in Figure 4. Figure 3 shows that the AV leaves some extra headway to tame the effect of the stop and go waves. Such waves continue to reappear but the control is able to dissipate partly the effects. The velocity profiles, Figure 4, show a strong oscillation reduction for the AV in the time interval [START_REF] Zhou | A recurrent neural network based microscopic car following model to predict traffic oscillation[END_REF]100], then oscillations affect the AV and finally are more under control after time t = 160.

Next we show that stop and go waves can be dampened using the control laws ( 13) and [START_REF] Gomes | Optimal freeway ramp metering using the asymmetric cell transmission model[END_REF]. Since the autonomous vehicle used in the experiment runs with velocity control, we will slightly modify our controls ( 13) and ( 15) as follows; let uj ∈ IR+ defined by uj = u(tj) with tj = j * ∆t, j ∈ IN and ∆t > 0. First, combining ( 13) with ( 7), we have vn = -(vn -v) with v ∈ IR+. Thus, the target velocity is

v target 1,j+1 = (uj -v) exp(-∆t) + v. ( 16
)
Since vn-1(s) = vn-1(t) := vn-1,j for every s ∈ [tj, tj + ∆t), combining ( 15) with ( 7), the target velocity is

v target 2,j+1 = uj - vn-1,j + v 2 exp(-∆t) + vn-1,j + v 2 . ( 17
)
From now on, v is a function depending on time and at every time t > 0, v(t) is constructed as the minimum of the temporal average speeds of the lead vehicle and the autonomous vehicle over (0, t). More precisely, the target velocity is vj = min Figure 5: Vehicles trajectories when using the control ( 19) with ( 16).

j-1 i=1 v lead i j -1 , j-1 i=1 ui j -1 (18)
Figure 6: Vehicles trajectories when using the control [START_REF] Hegyi | Optimal coordination of variable speed limits to suppress shock waves[END_REF] with [START_REF] Han | Resolving freeway jam waves by discrete first-order model-based predictive control of variable speed limits[END_REF].

The theory described in Section 4.2 is still useful because, for j large enough, vj ≈ c with c > 0 a constant. Adding the rule [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF] described in Section 4.1, the commanded velocities for ( 16) and ( 17) are updated via uj+1 = βj(αjv target k,j

+ (1 -αj)v lead j ) + (1 -βj)uj , k ∈ {1, 2}. ( 19
)
where v target 1,j and v target 2,j are defined in ( 16) and ( 17) respectively replacing v by vj defined in [START_REF] Hegyi | Model predictive control for optimal coordination of ramp metering and variable speed limits[END_REF]. Above, αj and βj are constructed as explained in Section 4.1. In Figure 5 and in Figure 6, we use the commanded velocity defined in [START_REF] Hegyi | Optimal coordination of variable speed limits to suppress shock waves[END_REF] with v target 1,j and v target 2,j respectively. Since the autonomous vehicle drives according to [START_REF] Hegyi | Model predictive control for optimal coordination of ramp metering and variable speed limits[END_REF], a gap is created when the lead car is affected by a stop and go wave.The difference between these controls and the P ID control, defined in Section 4.1, is that the autonomous vehicle never needs to catch up the lead vehicle and therefore it doesn't create another stop and go wave for the controls defined in [START_REF] Hegyi | Optimal coordination of variable speed limits to suppress shock waves[END_REF]. More precisely, the term min(max( ∆x-7m 23m , 0), 1) in ( 10) is not needed anymore. Moreover, in [START_REF] Hegyi | Model predictive control for optimal coordination of ramp metering and variable speed limits[END_REF], we don't use that the length of the ring is equal to L.

The effect of the Lyapunov controllers can be seen in the velocity profiles in Figures 7,8 which shows that the traffic is smoother when the control is active. More precisely, in Figure 7, we notice that the AV tends to have very small oscillations when the control is active. On the other side, the not controlled vehicles tend to keep some velocity oscillations, however apparently less than the PID controller. The interesting fact about this control is that the system, once introduced the saturation, tends to find local equilibrium and with some initial data might only dampen the stop and go waves without dissipating them. Figure 8 show a behavior pretty similar to that of the PID controller, with even stronger oscillations for the AV.

Fuel consumption

Fuel consumption is estimated for the simulation results using the VT-Micro fuel consumption model [START_REF] Rakha | Comparison of mobile5a, mobile6, vt-micro, and cmem models for estimating hot-stabilized light-duty gasoline vehicle emissions[END_REF][START_REF] Rakha | Development of VT-Micro model for estimating hot stabilized light duty vehicle and truck emissions[END_REF]. The VT-Micro model uses polynomial regression on vehicle fuel consumption data collected at Oak Ridge National Lab to predict fuel consumption based on vehicle speed and acceleration. In simulation, per-vehicle fuel consumption averaged over an ensemble of 10 simulation runs when all traffic is under human control is 21.93 /100 km. When the PID controller ( 11) is implemented on the AV in simulation and the remaining 21 vehicles are under human control, the average per-vehicle fuel consumption is reduced by 16.3% to 18.36 /100 km. Similarly, when the controller in ( 19) with ( 16) is used on the AV, the average per-vehicle fuel consumption is reduced by 17.8% to 18.03 /100 km, while when the controller in ( 19) with ( 17) is used on the AV, the average per-vehicle fuel consumption is reduced by 17.3% to 18.13 /100 km. The fuel consumption means over the 10 simulation runs conducted and corresponding standard deviations are presented in Table 1. As described in Section 4.3, the Lyapunov controllers might not always dissipate the stop-and-go waves which results in a higher standard deviation for the mean fuel consumption in Table 1 for these controls.

5 Experimental results

Experimental design

A series of experiments were conducted to validate the development of stop and go traffic waves and to demonstrate the dissipation of these waves using the designed controller. We first describe the experiment setup and its benefits and limitations. Then we present the experimental results of a low penetration rate of autonomous vehicles actively dissipating stop and go waves. The experiments are conducted using a similar setup to the seminal works of Sugiyama, et al. [START_REF] Sugiyama | Traffic jams without bottlenecks-experimental evidence for the physical mechanism of the formation of a jam[END_REF] and Tadaki, et al. [START_REF] S.-I. Tadaki | Phase transition in traffic jam experiment on a circuit[END_REF], which was able to isolate human driving behavior without considering other factors. These experiments showed conclusively that human driving behavior alone is sufficient to trigger stop and go waves. In this experiment, we use the same setup because it has been shown to reliably produce the kinds of stop and go waves that the controllers in this paper are designed to dampen.

While the underlying experimental setup in the work presented is similar to the setup used by Sugiyama, et al. [START_REF] Sugiyama | Traffic jams without bottlenecks-experimental evidence for the physical mechanism of the formation of a jam[END_REF], minor modifications were made to accommodate for the larger US vehicles. Just as in [START_REF] Sugiyama | Traffic jams without bottlenecks-experimental evidence for the physical mechanism of the formation of a jam[END_REF], the experiment is conducted with a total of 22 vehicles on the track. However, due to the substantially larger US vehicles, the track length was increased to 260 m around. While all vehicles have human-drivers in them, one of the vehicles, the University of Arizona Cognitive Autonomous Test Vehicle (CAT Vehicle) is an autonomous-capable vehicle and can be switched from being human-piloted to autonomous during the experiment.

The experiment begins with all vehicles evenly spaced and at rest. When the drivers are given a signal, they begin to drive, and human-piloted traffic conditions are observed. After some time the Vehicle trajectories of each vehicle were collected using a VSN Mobile V360 panoramic video camera placed at the center of the track. The video footage recorded during the experiment was processed using image processing algorithms. More details on the image processing algorithms used can be found in the article by Wu, et al. [START_REF] Wu | Tracking vehicle trajectories and fuel rates in oscillatory traffic[END_REF]. Additionally, vehicle performance data such as fuel consumption was recorded during the experiment using OBDLink MX onboard diagnostics (OBD-II) data loggers.

The ring-road experimental design is selected because it has been shown to produce traffic instabilities that are similar to those observed in real highway traffic [START_REF] Sugiyama | Traffic jams without bottlenecks-experimental evidence for the physical mechanism of the formation of a jam[END_REF]. Furthermore, using a closed-circuit experimental test track allowed us to reproduce "infinite" traffic (where each vehicle has a vehicle in front of it and a vehicle behind it at all times) with a finite number of vehicles. However, there are some limitations that arise due to this experimental design. Since the track is only a single lane of traffic, this experimental setup cannot be used to assess the robustness of the designed algorithm to overtaking and merging.

Experimental results

An experiment was conducted to test the PID velocity controller in a field test. The experiment was conducted on a large flat parking lot in Tucson, AZ. The human-piloted vehicles used during this experiment were rented from the University of Arizona Motor Pool, while the AV that was used was the CAT Vehicle. An overview of the test track is seen in Figure 9 where all 22 vehicles are on the track following the experiment. The AV used to implement the traffic controller, the University of Arizona CAT Vehicle, is seen in Figure 10 being parked at the start position for the experiment.

The experiment is started with all vehicles at rest on the track and under human control. The velocity profile and vehicle trajectories for all vehicles in the experiment is presented in Figure 11 and Figure 12, respectively. Here, the CAT Vehicles speed and trajectory are plotted in red, while the speed and trajectory for the remaining vehicles are plotted in grey. After 161 seconds the small oscillations grow, and a noticeable stop and go wave develops. In the presence of this wave, vehicles fluctuate in speed between 0 m/s and 14 m/s as seen in Figure 11. This same stop and go wave is also seen in the vehicle trajectories in Figure 12. After 218 seconds the controller on the CAT Vehicle is activated and the traffic is under AV control. This control action is maintained until the experiment is ended after 413 seconds.

As seen in the vehicle trajectories in Figure 12, after roughly 320 seconds, one of the human-piloted vehicles introduces a small traffic wave. This wave propagates upstream until it encounters the CAT Vehicle, which is able to partially dampen, but not fully remove the wave. On the second pass around the track, the CAT Vehicle is able to fully eliminate the wave. This demonstrates that a single AV may not be sufficient to fully eliminate stop and go waves, but a low penetration rate of AVs may be sufficient to eliminate, or at least substantially dampen stop and go waves. The effect that the traffic controller implemented on the single AV has on the overall traffic flow (all 21 human-piloted vehicles in addition to the CAT Vehicle) is apparent in the velocity profile in Figure 11 and vehicle trajectories 12 where substantially smoother traffic is seen when the traffic is under the control of the AV. Specifically, the average velocity standard deviation between the time period where waves are present and the control period is reduced by 54.7% from 3.85 m/s to 1.74 m/s while not substantially effecting throughput (1755 veh/hr when waves are present and 1711 veh/hr when the AV is actively dissipating the stop and go waves, or a reduction of 2.5%). Furthermore, since each vehicle in the flow was instrumented with an OBD-II scanner, it was possible to measure the instantaneous fuel consumption of each vehicle during the experiment. The average fuel consumption over the entire fleet when waves are present was 26.3 /100 km, while the average fuel consumption of all vehicles when the AV is actively dissipating stop and go waves is reduced by 20.7% to 26.3 /100 km. It is important to note that this reduction in fuel consumption is not only realized on the CAT Vehicle, but is the averaged quantity across all vehicles in the experiment.

This experimental result validates the theoretical and simulation-based results and demonstrates that even at a low penetration rate, AVs are capable of substantially improving traffic flow by reducing velocity standard deviation. This leads to smoother traffic and reduces fuel consumption of not only the AV, but all the vehicles in the traffic flow.

Conclusion

In this work we establish controllability results for two well known microscopic traffic flow models in the setting in which an autonomous vehicle is able to be controlled with the aim of dissipating stop and go traffic waves. Based on the controllability results, three control algorithms are developed and tested in simulation, indicating that fuel consumption reductions of up to 20% may be achieved when the flow is stabilized by the autonomous vehicle. In a first of its kind field test, we further establish for one of the proposed control algorithms that traffic wave dissipation is possible with real autonomous and human piloted vehicles, and leads to a substantial reduction in fuel consumption compared to when waves are present.

While our results show the feasibility of control of traffic via AV, especially for fuel consumption reduction, there are some limitations to our study. For instance multilane traffic was not considered and, more generally, we did not include additional challenges coming from the impact of combined complex phenomena of multilane dynamics, merges, ramps, and non-FIFO assumptions.

Looking forward, we note that the level of difficulty to prove controllability results of the combined optimal velocity follow the leader model for n > 9 suggests alternate modeling scales, e.g., in the micromacro direction [START_REF] Delle Monache | Scalar conservation laws with moving constraints arising in traffic flow modeling: an existence result[END_REF][START_REF] Chalons | A conservative scheme for non-classical solutions to a strongly coupled pde-ode problem[END_REF][START_REF] Lattanzio | Moving bottlenecks in car traffic flow: A PDE-ODE coupled model[END_REF][START_REF] Piacentini | Traffic control via moving bottleneck of coordinated vehicles[END_REF], might be promising.

                      . ( 21
)
For n ∈ IN define the matrices (Dn-1, Fn-1) ∈ M2n-1(IR) 2 by:

D n-1 = -γ γ (0) -γ γ (0) -γ               , F n-1 = -(α + b) α (0) -(α + b) α (0) -(α + b)               . ( 22
)
with γ = bV (y * + d) and α = a (y * +d) 2 . Moreover, denote by 0n-1 ∈ Mn-1(IR) and In-1 ∈ Mn-1(IR) the zero matrix and the identity matrix respectively. Let k 1. From ( 21), there exist (A i,k ) i∈{1,••• ,4} ∈ Mn-1(IR) and (C i,k ) i∈{1,2} ∈ Mn-1,1(IR) such that,

A k =   A 1,k A 2,k C 1,k A 3,k A 4,k C 2,
Combining ( 23) with [START_REF] Jurdjevic | Controllability and stability[END_REF] and using that, for every k 1, A k B(2n -1) = 0, we conclude that, for every k 2,

A k+1 B = DA k-1 B + FA k B, (25)
Moreover, we have

A 2 B = X + F2n-1AB with AB =         0n-2,1 1 0n-3,1 α -α 0         and X =         0n-2,1 b 0n-3,1 γ -γ 0         . (27)

Figure 1 :

 1 Figure 1: Sketch of n vehicles on a trajectory along a ring of length L.

Figure 2 :

 2 Figure 2: Vehicle trajectories with no control. Note the red trajectory follows a single vehicle which acts identically to all grey trajectories.

Figure 3 :

 3 Figure 3: Vehicle trajectories of uncontrolled vehicles (grey) when using the PID control (11) applied only to a single vehicle (red).

Figure 4 :

 4 Figure 4: Velocity profiles of uncontrolled vehicles (grey) when applying the PID control to a single vehicle (red) (11).

Figure 7 :

 7 Figure 7: Simulated velocity profiles of human-piloted vehicles (grey) and autonomous vehicle (red) when using the control (19) with (16).

Figure 8 :

 8 Figure 8: Simulated velocity profiles of human-piloted vehicles (grey) and autonomous vehicle (red) when using the control (19) with (17).

Figure 9 :

 9 Figure 9: Vehicles on test track during the experiment.

Figure 10 :

 10 Figure 10: The CAT Vehicle, the autonomous vehicle used during the experiment.

Figure 11 :

 11 Figure 11: Velocity profile of all vehicles in the experiment. The grey velocity profiles are the human-piloted vehicles and the red profile is the CAT Vehicle. The blue vertical lines mark key times during the experiment: the start, at what point waves are clearly visible, when control of the CAT Vehicle is activated, and the end of the experiment.

Figure 12 :

 12 Figure 12: Trajectories of all vehicles in the experiment. The grey trajectories are the human-piloted vehicles and the red trajectory is the CAT Vehicle. The blue vertical lines mark key times during the experiment: the start, at what point waves are clearly visible, when control of the CAT Vehicle is activated, and the end of the experiment.

  . . . . . . . . . . . . . . . . 0

A 4 ,A 4 ,

 44 k+1 = Dn-1A 2,k + Fn-1A 4,k , A 2,k+1 = A 4,k , C 2,k+1 = Dn-1C 1,k + Fn-1C 2,k , C 1,k+1 = C 2,k . k+1 = Dn-1A 4,k-1 + Fn-1A 4,k , A 2,k+1 = Dn-1A 2,k-1 + Fn-1A 2,k , C 2,k+1 = Dn-1C 2,k-1 + Fn-1C 2,k , C 2,k = Dn-1C 1,k-1 + Fn-1C 1,k .

Table 1 :

 1 Simulation fuel consumption estimates over an ensemble of 10 simulation runs using the VT-Micro fuel consumption model.

	Control	Mean fuel consumption ( /100 km)	Standard deviation ( /100 km)
	No control	21.93	0.94
	PID controller (11)	18.36	0.35
	(19) with (16)	18.03	0.46
	(19) with (17)	18.13	0.65

Established for n <= 9 vehicles and conjectured for n > 9.

Acknowledgement

This research was supported by the Inria associated team "ModEling autonoMous vEhicles iN Traffic flOw " (MEMENTO). This material is based upon work supported by the National Science Foundation under Grant No. CNS-1446715 (Piccoli), CNS-1446690 (Seibold), CNS-1446435 (Sprinkle), and CNS-1446702 (Work). The authors thank Hannah Pohlmann for her help in designing and tuning of the Follow-the-Leader-Bando microscopic model.

Appendix

The linearization of system [START_REF] Coron | Control and nonlinearity[END_REF] at the equilibrium point (Y * , 0) ∈ IR 2n is described by

Proof of Theorem 3. If b = 0 then γ = 0. From [START_REF] Khodayari | A modified car-following model based on a neural network model of the human driver effects[END_REF] and [START_REF] Konishi | Decentralized delayed-feedback control of an optimal velocity traffic model[END_REF], for every k 2,

Thus, the Kalman controllability matrix (6) satisfies

By Cayley-Hamilton Theorem, there exists (α0,

Using the expression of AB given in [START_REF] Konishi | Decentralized delayed-feedback control of an optimal velocity traffic model[END_REF] and the equality

whence the conclusion.

Proof of Theorem 4. Let's prove by induction that, for every k 1, there exist (λi

.

Thus, (P k ) holds for k = 1. Assuming that (P k ) holds for k = p. From (25), we have

Using ( 28) and (P k ) for k = p, we conclude that (P k ) holds for k = p + 1.

The equality (P k ) for k = n -2 gives

Since AB and A 2 B are linearly independent and

Thus, the linearization of system (8) at the equilibrium point (Y * , 0) is controllable. Using [8, Theorem 3.8], Theorem 4 is proved.

Proof of Theorem 5. Using using the symbolic mathematics software Maple [https://www.maplesoft. com/products/Maple/], we establish that for every 3 n 9,

.

Thus, for every α = γ b , γ = 0, b = 0, the linearization of system (8) at the equilibrium point (Y * , 0) is controllable, whence the conclusion of Theorem 5 by using [START_REF] Coron | Control and nonlinearity[END_REF]Theorem 3.8].


Bertha Mayela

Toledo Acosta

Xavier Heiligenstein

Grégoire Malandain

Patrick Bouthemy

INTENSITY-BASED MATCHING AND REGISTRATION FOR 3D CORRELATIVE MICROSCOPY WITH LARGE DISCREPANCIES

Keywords: 3D microscopy images, 3D CLEM, multimodal patch matching, 3D registration

Correlative microscopy, especially light and electron microscopy (CLEM), enables the study of cells and subcellular elements in complementary ways, provided a reliable registration between images is efficiently achievable. We propose a general automatic registration method. Due to large discrepancies in appearance, field-of-view, resolution and position, a pre-alignment stage is required before any 3D fine registration stage. We define an intensity-based method for both stages, which leverages a common representation of the two involved image modalities. We report experimental results on different real datasets of 3D correlative microscopy, demonstrating time efficiency and overlay accuracy.

INTRODUCTION

In recent years, correlative microscopy has become an important and powerful tool in the bioimaging field. It combines two (or more) microscopy modalities, allowing the study of biological specimens in complementary ways and at different resolutions. Since current tools used by biologists require manual intervention, one key challenge in this area is to make the image overlay between modalities as automatic as possible for an easier and more efficient workflow. It will also favour dissemination of correlative microscopy and enable handling of large amount of images. This is particularly crucial for 3D acquisition. Our goal is precisely to propose an automated registration method for 3D correlative microscopy.

3D multimodal image registration has a long story in medical imaging [START_REF] Maintz | A survey of medical image registration[END_REF]. Briefly speaking, there are two main approaches for registration [START_REF] Zitová | Image registration methods: a survey[END_REF]. The first one is the geometric approach which first extracts specific points in each image, pairs (or matches) them, and computes the 3D registration transformation from the matched pairs of points. It can handle large transformations. The second approach is the intensity-based one, which iteratively computes the 3D transformation by minimizing a similarity criterion for instance based on mutual information. It implicitly assumes that an initial positioning is given, not too far from the optimal solution. However, 3D correlative microscopy raises specific issues, which makes both approaches not straightforwardly applicable for an automatic workflow as explained just below.

The most frequent combination of correlative microscopy is between light microscopy (LM) and electron microscopy (EM), referred to as correlative light and electron microscopy (CLEM). Fluorescence LM allows identification and tracking of labeled biomolecules, and provides information on their dynamics and interactions. However, the diffraction limits the resolution, and fluorescence LM cannot visualize unlabeled cellular structures. EM displays much higher resolution, and it can delineate cellular ultrastructure but at the cost of fixing the sample. Since there are several types of LM and EM modalities, any combination results in a specific CLEM method, also dependent on the biological problems under study. Other imaging techniques may accompany CLEM, such as microCT, or X-rays [START_REF] Schaad | Correlative imaging of the murine hind limb vasculature and muscle tissue by MicroCT and light microscopy[END_REF][START_REF] Handschuh | A correlative approach for combining microCT, light and transmission electron microscopy in a single 3D scenario[END_REF][START_REF] Karreman | Chapter 13 -find your way with X-Ray: Using microCT to correlate in vivo imaging with 3D electron microscopy[END_REF].

The technological and technical differences between EM and LM images are in many aspects: field of view, pixel resolution, image size, content, appearance, raising important and specific challenges for automated registration in 3D correlative microscopy. Furthermore, large gaps occur between the initial locations and orientations of the objects of interest in the respective LM and EM stacks in the data acquisition workflow. Then, geometrical and appearance discrepancies between LM and EM stacks hamper a direct computation of the 3D transformation mapping one image stack to the other.

The geometric registration approach is the one usually adopted in the semi-automatic 3D correlative microscopy workflow, but the selection and pairing of corresponding points are typically manually performed. In contrast, existing automatic segmentation methods cannot deliver repeatable point extraction in all LM and EM images. Indeed, these methods may not extract a sufficient amount of reliable point pairs in the two image modalities, in particular due to the lack of structure and texture in fluorescence images. The intensity-based approach for direct 3D registration computation is usually employed in medical imaging, since it can benefit from close initialisation and isotropic volumes. This does not stand for correlative microscopy which then requires a pre-alignment stage [START_REF] Kopek | Diverse protocols for correlative superresolution fluorescence imaging and electron microscopy of chemically fixed samples[END_REF].

The paper is organized as follows. In Section 2, we briefly survey related work in CLEM. Section 3 introduces the datasets and the corresponding image properties. In Section 4, we present the two stages of our method. Section 5 reports experimental results, and provides a comparison with manual interaction using eC-CLEM [START_REF] Paul-Gilloteaux | eC-CLEM: flexible multidimensional registration software for correlative microscopies[END_REF]. Section 6 contains concluding remarks.

RELATED WORK IN CLEM

In current CLEM setups, registration is performed by a first global pre-alignment followed by a refined alignment, and usually involved an important user interaction. The manual pre-alignment or matching of the EM and LM stacks is commonly performed by biologists, either during the acquisition step and/or during the post-acquisition registration step. During acquisition, the pre-alignment relies on vi-sually identifying predefined marks on the sample holder, or exploiting fiducial markers introduced in the sample [START_REF] Peddie | Correlative and integrated light and electron microscopy of in-resin GFP fluorescence, used to localise diacylglycerol in mammalian cells[END_REF][START_REF] Franks | Correlative Fluorescence and Electron Microscopy in 3D-Scanning Electron Microscope Perspective[END_REF]. Likewise, during the registration stage, the biologist has to select pairs of points in the two images, whatever they are (fiducial markers, distinctive features on the sample) [START_REF] Beckwith | Seeing a Mycobacterium-Infected cell in nanoscale 3D: Correlative imaging by light microscopy and FIB/SEM tomography[END_REF][START_REF] Cortese | Chapter 6 -3D HDO-CLEM: Cellular compartment analysis by correlative Light-Electron microscopy on cryosection[END_REF]. The landmark pairs are then exploited to compute the 3D geometric transformation between EM and LM stacks [START_REF] Fukuda | Coordinate transformation based cryo-correlative methods for electron tomography and focused ion beam milling[END_REF][START_REF] Russell | 3D correlative light and electron microscopy of cultured cells using serial blockface scanning electron microscopy[END_REF][START_REF] Bohórquez | Correlative confocal and 3D electron microscopy of a specific sensory cell[END_REF]. The matching step is a demanding and tedious task, especially if a large amount of data must be processed [START_REF] Maco | Chapter 16 -correlative in vivo 2-Photon imaging and focused ion beam scanning electron microscopy: 3D analysis of neuronal ultrastructure[END_REF]. Such semi-automatic framework is adopted in [START_REF] Handschuh | A correlative approach for combining microCT, light and transmission electron microscopy in a single 3D scenario[END_REF] for registering 3D microCT and LM images, in the recently released eC-CLEM software for 3D correlative microscopy [START_REF] Paul-Gilloteaux | eC-CLEM: flexible multidimensional registration software for correlative microscopies[END_REF], or in [START_REF] Murphy | Correlative 3D imaging of whole mammalian cells with light and electron microscopy[END_REF] to overlay Scanning EM (SEM) and confocal LM volumes.

In above works, the refined registration is addressed by pairing corresponding points, and the user has to iteratively re-apply this procedure (picking points, computing the transformation) until he or she found the computed transformation satisfactory. A few works propose an automated refined registration (the pre-alignment being still manual). In [START_REF] Murphy | Correlative 3D imaging of whole mammalian cells with light and electron microscopy[END_REF], the cell membrane is segmented in the SEM stack and modified to a resolution similar to LM one. Then, EM and LM volumes are registered with an affine transformation and normalized correlation. [START_REF] López | Chapter 8 -a fully integrated, three-dimensional fluorescence to electron microscopy correlative workflow[END_REF] proposes an automated refined registration based on mutual information and rigid transformation. In [START_REF] Paul-Gilloteaux | eC-CLEM: flexible multidimensional registration software for correlative microscopies[END_REF], an extension to an automatic registration was also proposed, provided spots can be correctly detected in both images, but it was only tested on 2D real CLEM.

Our contribution is to propose automated procedures for both the pre-and the refined alignment for 3D correlative microscopy. Our method is intensity-based for both stages, and leverages a common representation of EM and LM images.

3D CLEM FRAMEWORK

To motivate our contribution, we first present the two 3D correlative microscopy sets we will deal with, and we give information on the acquired images. The first set is composed of a serial block face SBFSEM image stack and a confocal image stack. HeLa cells were imaged for LM in a 710 LSM Zeiss inverted confocal microscope; afterwards, they were fixed, stained and embedded for SBF-SEM imaging and acquired using a Sigma VP scanning electron microscope (Carl Zeiss). The LM stacks has two channels: transmitted light and one fluorescent channel, with volume size of 1024×1024×16 and pixel size of 0.22µm in X and Y , and 0.6µm in Z, knowing that the XY -space corresponds to the lateral dimensions of the stack and the Z-axis to the vertical dimension of the stack. The size of the EM stack is 666 × 903 × 247 and its pixel size is 0.14µm in X and Y , and 50nm in Z, which is rather lowmagnification EM. A 3D rendering is displayed in In both datasets, the 3D stacks are anisotropic and have significant differences in terms of volume size, pixel resolution, content, appearance, and location of structures of interest between modalities. Nevertheless, in these examples, it turns out that the discrepancy between the two volumes in terms of shift along the Z-axis and of rotation in the XY plane are not that significant, but other discrepancies remain.

OUR AUTOMATED OVERLAY METHOD

Intensity-based XY pre-alignment

To address the large gap between the initial positions of EM and LM stacks, we have defined an intensity-based prealignment method. We restrict to the shift in the XY -space, since it is the larger component of misalignment (which is often the case by the way). It is tractable for an automatic pre-alignment at that stage. To compute this initial XY shift, we accordingly reduce the dataset to two dimensions. Projection along the Z-axis of the LM stack, using the maximum intensity projection (MIP), yields a good distribution of the volume content. If the LM stack is of size m × n × l, LM-MIP is a 2D m × n image obtained by scanning LM intensities along the Z-axis and selecting the highest one for each (x, y) location of the image grid. The same cannot be applied to the EM stack, since the resulting MIP image will be saturated. Then, we simply select a few slices in the EM stack. CLEM overlay may be performed in both ways, EM to LM, or LM to EM. The overlay is focused on a region of interest (ROI) of the biological sample, which may involve a specific cell or a subcellular structure in a cell, and is specific to each biological project. The ROI is expected to be delineated either by the biologist (which will be in that case the only user interaction), or by an algorithm able to detect the content of interest.

Here we start with the EM stack. Given a 3D ROI delineated in the EM stack, we take the corresponding 2D ROIs in the selected EM slices, and we search for their corresponding ROI in the LM-MPI image. The patch dimensions in the LM-MPI image are inferred from the EM-ROI ones, while taking into account the ratio between pixel resolutions available in the metadata attached to the acquisition. To facilitate the search, we require a comparable appearance for the two images. To do so, we leverage the Laplacian of Gaussian (LoG) transform we introduced in [START_REF] Toledo-Acosta | A common image representation and a patch-based search for correlative lightelectron-microscopy (CLEM) registration[END_REF] for 2D CLEM registration. The LoG transform maps a 2D EM image and a 2D LM image to closer representation, decreasing the big difference in appearance between the images as illustrated in Fig. 4. The LoG transform enhances high spatial frequencies while suppressing local linear intensity variations. The LoG transform of an image I is expressed as:

Lσ(I) = ∆(Gσ * I), (1)
where the Gaussian standard deviation σ acts as an adaptable scale linked to the size of the objects of interest in I. It is automatically selected by using the method described in [START_REF] Basset | Adaptive spot detection with optimal scale selection in fluorescence microscopy images[END_REF] which exploits the Lindeberg's scale-space approach. We select the scale with the maximum number of blobs, which are local minima in the constructed scale-space domain. By adapting the LoG scale to the EM and LM image content, we lessen the difference in content scale of LM and EM images. Once a common representation of EM and LM images is created, we have to match the 2D LoG-EM-ROI with its corresponding ROI in the 2D LoG-LM-MIP image. We handle several 2D LoG-EM-ROIs, one by EM slice, to ensure robustness by making the matching not dependent on one particular 2D EM slice. Then, a 2D patch-based exhaustive search is carried out over the LoG-LM-MIP. Due to the differences in orientation between the LM and EM images, it is necessary to adopt a similarity measure that is invariant to rotation and scaling. Therefore, we take histogram-based descriptors for the matching process. We have used two types of histograms: the LoG-values histogram and the Local Directional Pattern (LDP) [START_REF] Jabid | Facial expression recognition using Local Directional Pattern (LDP)[END_REF]. We have also implemented two histogram distances as similarity measures: the cosine similarity and the histogram intersection distance to evaluate the similarity of the LoG-EM-ROI with each tested LoG-MIP-LM patch.

If we use n EM slices taken from the EM stack, we have n candidate patches for the ROI location in the 2D LM-MIP image as illustrated in Fig. 3. The final ROI center location in the LM-MIP image is selected using a robust criterion since outliers may occur, i.e., wrong matching. It is given by the weighted average of the center positions of the inlier patch candidates. The weights are defined as an inverse of the distance between the center of the patch and the median of the centers positions of all candidate patches. Beforehand, outliers are discarded according to the median absolute deviation of their center position to the median of the centers of all selected patches. The difference between the center location of the Fig. 3. ROIs extracted from multiple LoG-EM slices to be located in the LoG-LM-MIP image selected ROI in the LoG-LM-MIP image and the center of any 2D EM-ROI (since the 2D EM-ROIs are vertically aligned in the 3D EM-ROI stack) yields the 2D shift in the XY -space. Then, the full LM stack is shifted by applying the estimated 2D translation.

To summarize, our pre-alignment method comprises four main steps: 1) Maximum intensity projection of the LM stack to a 2D image; 2) LoG representation of 2D-LM-MPI image and 2D EM slices; 3) Matching of the ROI of every LoG-EM slice in LoG-LM-MIP image; 4) Pre-alignment of EM and LM stacks by a XY -shift.

3D ROI-based affine registration

Once the pre-alignment is achieved, we can estimate the 3D transformation which will allow to overlay the 3D LM stack onto the 3D EM stack. The 3D transformation is computed between the 3D LoG-EM-ROI and the pre-aligned 3D LoG-LM-ROI. The transformed LoG-LM stack is resampled to fit the same size in the Z-dimension as the LoG-EM stack using a bilinear interpolator after applying a low-pass filter. The 3D LoG-LM-ROI is recovered using the same coordinates of the 3D LoG-EM-ROI.

The LoG images are still different enough so that an intensitydifference similarity measure cannot provide satisfactory results. We resort to mutual information (MI) as a similarity measure, due to its well-known capability of handling multimodal images. First, a 3D rigid transformation (composed of a rotation and a translation is estimated. However, given the 3D CLEM workflow, it is not sufficient to account for the geometric relationship between the LM and EM stacks. The EM slices do not lie in the same plane as the LM slices, and distortions between stacks may occur during the acquisition. Consequently, we refine the registration with a 3D affine transform.

Our overall method does not involve any extraction of image features from the EM and LM stacks. It is merely intensity-based, which prevents from any possible segmentation errors and saves computation time. Furthermore, it does not require any parameter setting by the user. We use ITK libraries in C++ for the MI-based registration.

EXPERIMENTAL RESULTS

Experiment 1: EM and LM. To pre-align the EM and LM stacks, we follow the 2D matching process described in Section 4.1. A 2D projection of the LM stack is generated using MIP. EM slices 107-118 were selected to extract a 2D ROI, whose LoG representation (Fig. 4a) is matched within the LoG-LM-MIP image. Using [START_REF] Basset | Adaptive spot detection with optimal scale selection in fluorescence microscopy images[END_REF], we obtain an approximate scale for the LoG transform of EM and LM, respectively estimated as σEM = 4.3 and σLM = 2.5.

A corresponding patch in the LoG-LM-MIP is found for each LoG-EM-ROI by the patch-based search, with a 90% overlap between tested patches, using both cosine and histogram intersection distances. Then, we robustly compute the selected patch from the weighted average of all their centers (Fig. 4b). The final overlay of the registered EM and LM ROIs is shown in Figure 5. The accuracy can be appreciated by observing the fluorescence (in green) correctly surrounding the cell nucleus in the EM slice. The computed 3D rigid registration corresponds to a 3D translation of (-1.21, -0.5, 7.17), and a rotation of angles -7.52 • , -2.63 • and -1.46 • around the X, Y and Z axis respectively. By picking corresponding points in LM and EM stacks, we evaluated an average registration error of 6.5 pixels, moderately less accurate than with manual eC-CLEM.

To compute this registration, it takes 1.54 minutes for the prealignment and 13.5 seconds for the 3D registration. Comparatively, manual correlation by a trained biologist on a known dataset using eC-CLEM took 20 minutes, most time being dedicated to primary orientation. Once the initial rotation, translation and scaling achieved (the first four points required in eC-CLEM), the fine correlation can be estimated to 5 minutes per points-pair. On an unknown dataset, initial correlation can take more than one hour.

CONCLUSION

We have defined an original and efficient intensity-based method for automatically overlaying 3D image stacks of different modalities in correlative microscopy (CM). Our method can cope with large appearance, resolution and position discrepancies. It proceeds in two stages: pre-alignment and 3D fine registration. It does not involve any parameter setting, and then, is straightforward to use by biologists. Results on different combinations of 3D correlative microscopy show its genericity. Registration accuracy is satisfying for a fully automatic overlay of the two 3D image stacks. Future work will deal with still more general handling of the pre-alignment (fully 3D shift) and possibly non-rigid 3D registration, while continuing to validate the proposed method on other types of 3D CM datasets.

Fig. 1 .

 1 Fig. 1. Exp. 1: 3D CLEM The second set, shown in Fig.2, includes a stack of microCT and a two-photon fluorescence microscopy stack. The microCT stack of

Fig. 4 .

 4 Fig. 4. Exp. 1: a) slice 107 of LoG-EM stack, b) LoG-LM-MIP image with patch candidates overprinted in magenta and the selected one in light bue (to be better visualized by zooming in the pdf file)

Experiment 2 :

 2 MicroCT and LM. In this experiment, slices 15-25 are extracted from the microCT stack, and a 2D ROI is extracted from each of their LoG representations. Then, a MIP of the LM stack is computed. The estimated LoG scales are σµCT = 22.2 and σLM = 2.1. The full LM stack is pre-aligned by applying the shift given by the first stage. Then, we compute the 3D affine transformation. Computation time is 3.9 minutes for the pre-alignment stage and 15.4 seconds for the 3D registration. Visual evaluation of the overlay is given in Fig.5. It is normal that most of the elongated shapes are not overlaid, since they are present in only one of the two modalities.

Fig. 5 .

 5 Fig. 5. 3D overlay for experiments 1 (left) and 2 (right). The fluorescence LM is overprinted in green. Top: Overlay visualized for one given slice. Bottom: Overlay visualized on the full 3D stack.

This work has been partly supported by the National Council of Science and Technology of Mexico (CONACYT) under the Grant no. 246372/384571.