Alberto Seeger
email: alberto.seeger@univ-avignon.fr

Jos É Vicente-P Érez

ON CARDINALITY OF PARETO SPECTRA *

Keywords: Cone-constrained eigenvalue problem, Pareto eigenvalue, Linear complementarity system

How many Pareto eigenvalues are there in a matrix of a prescribed order? This note provides the best lower bound that it is known up to now for the maximal number of Pareto eigenvalues in a matrix of order n.

1. Introduction. Throughout this work, the space R n is equipped with the usual inner product •, • and componentwise ordering. So, the notation x ≥ 0 indicates that each component of x is nonnegative. The dimension of the zero vector 0 is clear from the context. Let M n denote the linear space of real matrices of order n. A real number λ is a Pareto eigenvalue of A ∈ M n if the linear complementarity system x ≥ 0, Axλx ≥ 0,

x, Axλx = 0 has a nonzero solution x ∈ R n . Such nonzero vector x is called a Pareto eigenvector of A. The set of all Pareto eigenvalues, denoted by Π(A), is called the Pareto spectrum of A. This terminology has been introduced in [START_REF] Seeger | Eigenvalue analysis of equilibrium processes defined by linear complementarity conditions[END_REF] and further employed in [START_REF] Adly | A nonsmooth algorithm for cone-constrained eigenvalue problems[END_REF][START_REF] Chen | The Pareto eigenvalue problem of some matrices (Chinese)[END_REF][START_REF] Hiriart-Urruty | A variational approach to copositive matrices[END_REF][START_REF] Da | Cone-constrained eigenvalue problems: theory and algorithms[END_REF]. Numerical methods for computing Pareto eigenvalues can be found, for example, in [START_REF] Adly | A nonsmooth algorithm for cone-constrained eigenvalue problems[END_REF][START_REF] Júdice | On the asymmetric eigenvalue complementarity problem[END_REF][START_REF] Júdice | The eigenvalue complementarity problem[END_REF][START_REF] Da | Cone-constrained eigenvalue problems: theory and algorithms[END_REF].

The concept of Pareto eigenvalue arises in various fields of applied mathematics.

Example 1.1. Consider an autonomous differential system of the form z(t) ≥ 0, Az(t)ż(t) ≥ 0, z(t), Az(t)ż(t) = 0.

A nonzero solution z(•) to the above system is obtained by setting z(t) = e λt x, where λ is a Pareto eigenvalue of A and x is an associated Pareto eigenvector.

The theory of Pareto spectra differs substantially from the classical spectral theory, in particular, when it comes to deal with cardinality issues. For instance, a matrix of order 10 has at most 10 eigenvalues, but the matrix can have more than 1500 Pareto eigenvalues! It remains a difficult open question to determine the exact value of

π n := max A∈Mn |Π(A)| ,
where |S| denotes the cardinality of a set S. The number π n is called the Pareto capacity of M n . It can be verified that π 1 = 1 and π 2 = 3, but π n for n ≥ 3 needs to be determined.

Example 1.2. A computation shows that

A =   1
3 9 -3 9 27 -9 27 81   has 9 Pareto eigenvalues. Hence, π 3 ≥ 9. One knows that a matrix of order 3 cannot have 11 Pareto eigenvalues or more (see [START_REF] Da | Cone-constrained eigenvalue problems: theory and algorithms[END_REF]), but it is not clear to us whether 10 Pareto eigenvalues is possible.

There is neither an exact nor an asymptotic formula that helps to determine π n . As shown in [START_REF] Da | Cone-constrained eigenvalue problems: theory and algorithms[END_REF], the sequence {π n } n≥1 is increasing and sandwiched as follows:

2 n -1 ≤ π n ≤ n(2 n-1 -1) + 1.
The above bounds are crude in general, and they are getting worse for larger values of n. The purpose of this note is to improve the lower bound on Pareto capacity π n .

2. Main result. For the sake of convenience we recall a useful lemma taken from [START_REF] Seeger | Eigenvalue analysis of equilibrium processes defined by linear complementarity conditions[END_REF] in the following. The notation A J = [a i,j ] i,j∈J denotes the principal submatrix of A which takes entries in the rows and columns of A indexed by J.

Lemma 2.1. The scalar λ ∈ R is a Pareto eigenvalue of A ∈ M n if and only if there exist a nonempty set J ⊂ {1, . . . , n} and a vector ξ ∈ R |J| such that

A J ξ = λξ, (2.1)
ξ j > 0 for all j ∈ J, (2.2) j∈J a i,j ξ j ≥ 0 for all i / ∈ J. (2.3)
Furthermore, a Pareto eigenvector x associated to λ is constructed by setting

x j = ξ j if j ∈ J, 0 if j / ∈ J.
Thus, in order to compute the whole Pareto spectrum of a matrix A ∈ M n , one must consider the 2 n -1 possible ways of selecting the index set J. For each choice of J one first solves the classical eigenvalue problem (2.1), and then checks whether or not the positivity constraint (2.2) and the binding constraint (2.3) are satisfied. The following theorem is our main result.

Theorem 2.2. For all n ≥ 2, 3(2 n-1 -1) ≤ π n . (2.4)
Proof. We construct a matrix of order n whose number of Pareto eigenvalues is equal to the left-hand side of (2.4). Consider the matrix

A =         s 2 s 3 s 4 s 5 • • • -s 3 s 4 s 5 s 6 • • • -s 4 s 5 s 6 s 7 • • • -s 5 s 6 s 7 s 8 • • • . . . . . . . . . . . . . . .         , (2.5) whose (i, j)-entry is a i,j = -s i+1 if i ≥ 2, j = 1, s i+j otherwise.
For a reason that will be clear in a moment, we force the parameter s ∈ R to satisfy the following two properties:

s > 1 + √ 2, (2.6) s is transcendental. (2.7)
In order to compute the Pareto spectrum of (2.5), we work out all the choices for the index set J = {j 1 , . . . , j q }. For convenience' sake, we split the job into several steps:

Step 1. We examine the case 1 / ∈ J. This case is easier because the condition (2.3) is automatically satisfied. By proceeding as in [7, Proposition 3] we get λ J := tr(A J ) = j∈J s 2j , ξ J := s j1 , . . . , s jq T as unique solution to (2.1)-(2.2). Strictly speaking, the eigenvector ξ J is unique only up to normalization. We deduce that λ J is a Pareto eigenvalue of A.

Step 2. We find the eigenvalues of A J when 1 ∈ J. If J = {1}, then A J has s 2 as unique eigenvalue. This subcase is of no interest because (2.2)-(2.3) becomes infeasible. Consider then 1 = j 1 < • • • < j q with q ≥ 2. We claim that A J has exactly two real nonzero eigenvalues. For proving this claim we expand the polynomial p J (λ) = det(A J -λE), where E denotes an identity matrix of appropriate size. The general theory of characteristic polynomials (cf. [4, Chapter 1]) asserts that

det(A -λE) = (-1) n λ n + n r=1 (-1) n-r c r (A)λ n-r with c r (A) =
|I|=r det(A I ). Applied to A J , the above determinantal expansion becomes

det(A J -λE) = (-1) q λ q + q r=1 (-1) q-r c r (A J )λ q-r with c r (A J ) = I⊂J, |I|=r det(A I ). (2.8)
While writing (2.8) we use the fact that (A J ) I = A I for I ⊂ J. Given the structure of (2.5), we see that det(A I ) = 0 if I contains at least two elements from {2, . . . , n}. Thus,

c 1 (A J ) = tr(A J ), c 2 (A J ) = j∈J\{1} det A {1,j} , c r (A J ) = 0 for all r ≥ 3 with det(A {1,j} ) = det a 1,1 a 1,j a j,1 a j,j = 2s 2(j+1) .
We now have p J (λ) = (-1) q λ q-2 Ψ J (λ) , where Ψ J (λ) := λ 2tr(A J )λ + w J , (2.9)

w J := j∈J\{1} 2s 2(j+1) = 2s 2 tr(A J ) -s 2 .
(2.10)

The root λ = 0 is thrown away because (2.1)-(2.2) becomes infeasible. The polynomial p J admits two other real roots, namely,

λ J ± = tr(A J ) ± √ ∆ J 2 .
The condition (2.6) yields tr(A J ) > 2(2 + √ 2)s 2 , which in turn implies the positivity of the discriminant

∆ J = tr(A J ) 2 -4w J = tr(A J ) 2 -8s 2 tr(A J ) + 8s 4 .
Step 3. We establish some useful bounds for λ J -and λ J + . We have

∆ J = s 2j2 + • • • + s 2jq + s 2 2 -8s 2 s 2j2 + • • • + s 2jq .
By expanding the square and rearranging terms, we get

∆ J = s 2j2 + • • • + s 2jq -s 2 2 -4s 2 s 2j2 + • • • + s 2jq ,
and therefore,

√ ∆ J < s 2j2 + • • • + s 2jq -s 2 = tr(A J ) -2s 2 .
This leads to the following chain of strict inequalities:

s 2 < λ J -< λ J + < tr(A J ) -s 2 . (2.11)
Step 4. We prove that λ J + and λ J -are Pareto eigenvalues of A. Consider the vectors ξ J -and ξ J + given by

ξ J ± = u J ± , s j2-n , s j3-n , . . . , s jq-n T with u J ± = tr(A J ) -s 2 λ J ± -s 2 s 1-n .
Due to (2.11), the term u J ± is well defined and positive. The remaining components of ξ J ± are also positive. Note that

(s 2 -λ J ± )u J ± + s 1+j2 s j2-n + • • • + s 1+jq s jq-n = -s 1-n tr(A J ) -s 2 + s 1-n s 2j2 + • • • + s 2jq = 0.
On the other hand, for all i ∈ {j 2 , . . . , j q }, we have

j∈J a i,j (ξ J ± ) j = -s 1+i u J ± + s i+j2 s j2-n + • • • + s i+jq s jq-n = -s 1+i u J ± + s i-n s 2j2 + • • • + s 2jq = -s 1+i tr(A J ) -s 2 λ J ± -s 2 s 1-n + s i-n tr(A J ) -s 2 = s i-n tr(A J ) -s 2 1 - s 2 λ J ± -s 2 = λ J ± s i-n = λ J ± (ξ J ± ) i .
The equality previous to the last one is due to the fact that

λ J ± 2 -tr(A J )λ J ± + 2s 2 tr(A J ) -s 2 = 0.
We have shown in this way that

A J ξ J ± = λ J ± ξ J ± . The formula j∈J a i,j (ξ J ± ) j = λ J ± s i-n (2.12)
applies also for i / ∈ J. Since the right-hand side of (2.12) is positive, the vector ξ J ± satisfies the inequality (2.3). This proves that λ J + and λ J -are Pareto eigenvalues of A.

Step 5. Finally, we determine the cardinality of Π(A). Both sets

J 1 = {J ⊂ {1, . . . , n} : 1 / ∈ J, J = ∅}, J 2 = {J ⊂ {1, . . . , n} : 1 ∈ J, |J| ≥ 2}
have the same cardinality, namely, 2 n-1 -1. Each J ∈ J 1 produces exactly one Pareto eigenvalue and each J ∈ J 2 produces exactly two Pareto eigenvalues. More precisely,

Π(A) = {λ J : J ∈ J 1 } ∪ {λ J -: J ∈ J 2 } ∪ {λ J + : J ∈ J 2 }. (2.13)
If all the elements on the right-hand side of (2.13) are distinct, then

|Π(A)| = 3(2 n-1 -1).
Let us examine whether there is a repetition in (2.13) or not. Take two distinct index sets I, J ∈ J 1 ∪ J 2 . There are three cases for consideration. The first case occurs when I and J are both in J 1 . The situation is then as in [7, Proposition 3] and we know that λ I = λ J . The second case occurs when I and J are both in J 2 . We claim that

{λ I -, λ I + } ∩ {λ J -, λ J + } = ∅. (2.14)
We prove that the polynomials Ψ I and Ψ J do not have a common root. Recall that Ψ J is given by (2.9)-(2.10). The polynomial Ψ I is defined, of course, in a similar way. Ab absurdo, let Ψ I (λ) = Ψ J (λ) = 0 for some λ ∈ R. Then λ 2tr(A I )λ + 2s 2 tr(A I )s 2 = 0, (2.15)

λ 2 -tr(A J )λ + 2s 2 tr(A J ) -s 2 = 0. (2.16)
The combination of (2.15) and (2.16) leads to λ tr(A J )tr(A I ) -2s 2 tr(A J )tr(A I ) = 0. But tr(A J )tr(A I ) = 0 because I = J. Hence, λ = 2s 2 . Substituting this value in (2.15) we get 4s 4 -2s 2 tr(A I ) + 2s 2 tr(A I )s 2 = 0. This leads to s = 0, which contradicts (2.6). The proof of (2.14) is complete. Consider now the third and last case: I ∈ J 1 and J ∈ J 2 . We claim that

λ I / ∈ {λ J -, λ J + }. (2.17)
The proof of this claim requires taking a look at the Pareto eigenvalues is a polynomial with integer coefficients. We shall prove that P I,J is non-constant. If P I,J were constant, then it should be identically zero because P I,J (s * ) = 0. Now, if P I,J is the zero function, then

λ I (s) = i∈I s 2i , λ J ± (s) = 1 2 j∈J s 2j ± 1 2      j∈J s 2j   2 -8s 2
λ I (s) ∈ {λ J -(s), λ J + (s)} for all s ∈ ]1+ √ 2, ∞[. By a continuity argument, either λ I (•) = λ J -(•) or λ I (•) = λ J + (•)
. But we have seen already that λ I (•) = λ J -(•) and λ I (•) = λ J + (•). Summarizing, s * is a root of a non-constant polynomial with integer coefficients. Hence, s * is algebraic. Therefore, the condition (2.17) holds true under the transcendentality (2.7).

Remark 2.3. If A s is the matrix given by (2.5), then |Π(A s )| = 3(2 n-1 -1) whenever s ∈ ]1 + √ 2, ∞[ is a transcendental number. This is enough for taking care of Theorem 2.2. The same proof shows that

S = s ∈]1 + √ 2, ∞[ : |Π(A s )| = 3(2 n-1 -1)
is a finite set containing only algebraic numbers. Numerical experiments suggest that S is empty, but this remains as a conjecture.

As one can see in Table 2.1, there is a large gap between the known lower and upper bounds for π n . Some work needs to be done to make this gap smaller, hoping to determine the exact value of π n .

  the parameter s. The term inside the square root is positive when s ranges on the open interval ]1+ √ 2, ∞[. On this interval the function λ J -(•) is infinitely often differentiable, but it is not a polynomial. Hence, λ I (•) and λ J -(•) are not the same function. For the same reason, the functions λ I (•) and λ J + (•) are different. Let us put ourselves in the worst situation and suppose that λ I (s * ) ∈ {λ J -(s * ), λ J + (s * )} (2.18) for a certain s * ∈ ]1 + √ 2, ∞[. The condition (2.18) says that Ψ J (λ I (s * )) = 0. Let us examine the term P I,J (s) = Ψ J (λ I (s)) as function of s ∈ ]1 + √ 2, ∞[. Clearly, P I,J (s) =

Table 2 . 1

 21 Lower and upper bounds for πn.

This author was supported by FPI Program of MICINN of Spain under grant BES-2006-14041.


Etienne Perot Valeo

Maximilian Jaritz
email: maximilian.jaritz@valeo.com

Raoul De Charette
email: raoul.de-charette@inria.fr

End-to-End Driving in a Realistic Racing Game with Deep Reinforcement Learning

Introduction

For autonomous driving the classical paradigm is to use a chain of perception, planning and control; but recent deep learning progress lets foresee an end-to-end alternative to map sensor inputs directly to low-level control of robots [START_REF] Levine | End-to-end training of deep visuomotor policies[END_REF]. End-to-end driving [START_REF] Mnih | Asynchronous methods for deep reinforcement learning[END_REF] was showcased in the car racing game TORCS using Reinforcement Learning but its physics and graphics lack realism.

We propose a method benefiting from latest asynchronous learning [START_REF] Mnih | Asynchronous methods for deep reinforcement learning[END_REF] to train an end-to-end agent in the context of a realistic car racing game -World Rally Championship 6 (WRC6). We do not rely on the in-game score and train solely on image and speed to learn the optimal action while reflecting real driving conditions. Our architecture was trained simultaneously on tracks with different graphics and road structure (cf. fig. 1 and3). Compared to previous use of TORCS [START_REF] Mnih | Asynchronous methods for deep reinforcement learning[END_REF][START_REF] Chen | Deepdriving: Learning affordance for direct perception in driving[END_REF][START_REF] Lau | Using Keras and Deep Deterministic Policy Gradient to play TORCS[END_REF], the environment exhibits more realistic physics (grip, drift), graphics (illuminations, animations, etc.), and a variety of environments (sharp turns, slopes, cliffs, snow, etc.). The proposed reward function converges faster than previous ones and offers some generalization capacity. Additionally, the driving style is more comparable to human driving.

Method

We used the asynchronous advantage actor-critic (A3C) [START_REF] Mnih | Asynchronous methods for deep reinforcement learning[END_REF] to train an end-to-end neural network. Every time-step, the algorithm receives the state of the game, acts (acceleration and steering), and gets a reward as supervision signal. This method optimizes driving policy using only RGB image as input (cf. fig. 1b) in order to maximize the cumulated reward. The choice of the A3C baseline is justified by its top performance and because it allows training without any need of experience replay for decorrelation. State encoder. Unlike other computer vision tasks a shallow CNN is sufficient as car racing relies mostly on road detection. Our CNN + LSTM architecture is similar to [START_REF] Kempka | A doom-based ai research platform for visual reinforcement learning[END_REF] though using a dense filtering (stride 1). It also uses max pooling to allow more translational invariance and takes advantage of speed and previous actions in the LSTM.

Reward shaping. To help the network to converge to the optimal set of solutions the reward shaping is crucial. The reason not to use the in-game score as reward is that the latter is too sparse to train the agents. In Mnih et al. [START_REF] Mnih | Asynchronous methods for deep reinforcement learning[END_REF] the reward R is computed with the angle difference θ between the car's heading and the road, and the speed v. Though efficient, it is limited as it does not prevent the car to slide along the guard rail since the latter follows the road angle.

Instead we chose to add the distance from the middle of the road d as a penalty, that is: R = v(cos θ -d).

Experiments

We ran the algorithm with 15 asynchronous agents. Each agent communicates via TCP with a WRC6 instance through a dedicated API specifically developed. It allows us to retrieve in-game info, compute the reward and send control back to the game. For computational reasons, costly graphics effects were removed and the horizontal field of view reduced. The game engine's clock runs at 30FPS and the physical engine is on hold as it waits for the next action.

Reward Comparison

To evaluate the proposed reward against the reward from [5], we trained a network with their architecture and plot the performance in fig. 2. Compared to [START_REF] Mnih | Asynchronous methods for deep reinforcement learning[END_REF] (blue curve), the proposed reward (green) converges faster (80% track completion after only 15 mega steps) while driving faster and safer. After 80 mega steps, the proposed reward drives at 88.0 km/h (i.e. +5.1km/h) and crashes 0.9 times per km (i.e. -5.3 crashes per km). The explanation is that with the previous reward the car tends to slide along the guard rail which slows it down and is more dangerous. Qualitatively also, the proposed reward leads to a smoother driving style. The agent had more difficulty to progress on the mountain track as it exhibits sharp curves and hairpin bends. cliffs, snow, etc. Physics also differ, especially road adherence. Guided back propagation is displayed in fig. 3 (i.e only positive inner gradients that lead to the chosen action). Despite the various scenes appearance, it learned to detect and use road edges/curvature as a strong control cue. This actually mimics classical approaches that also use lane markings for lateral controls. The training performance is also shown in fig. 4 and demonstrates that it progressed well in each track. After 190 mega steps the agent learned to drive in mountain, snow and coast tracks, to take some sharp turns and even hairpin bends. This is visible in fig. 1c that depicts crash along the tracks. We also tested the generalization of this training on unseen test tracks. It exhibits generalization capabilities as the bot was able to drive new tracks some An online video illustrates the over trained and new tracks: https://youtu.be/e9jk-lBWFlw

Figure 1 :

 1 Figure 1: The racing environment for our end-to-end driving architecture. (a) The full render in WRC6 game. (b)The 84x84 network input with guided back-propagation (blue overlay). Note, the narrow field of view and removal of indicators (turns signs, scores, etc.). End-to-end performance after training is displayed over tracks in (c).

Figure 2 :

 2 Figure 2: Proposed reward (green) versus reward from [5] (blue). For each reward the rolling mean (dark) and standard deviation (light) are shown (rolling uses 200 steps).

Figure 3 :

 3 Figure 3: Guided back propagation (blue highlights) of CNN+LSTM architecture after 190 mega steps. Despite the various scene and road appearances the network learned to detect road edges and to rely on the later for control.

  Performance and generalization A deeper network was trained for 190 mega steps, with three very different tracks (5 instances of each). The tracks contain sharp curves,

Figure 4 :

 4 Figure 4: Performance for the challenging training tracks.The agent had more difficulty to progress on the mountain track as it exhibits sharp curves and hairpin bends.


Sule Kahraman

Raoul De Charette

Influence of Fog on Computer Vision Algorithms

Keywords: Fog, Bad weather, Vision in bad weather, Computer Vision, Computer Graphics

This technical report describes a new preliminary approach to simulate fog in images using accurate physical and photometric models to study the influence of small particles on computer vision algorithms.

I. INTRODUCTION

The Earth's atmosphere is significantly non transparent in bad weather conditions. The small particles in the atmosphere create a scattering medium like fog that attenuates the light passing through and result in failure of computer vision algorithms such as feature detection, object recognition, tracking.

Particles with size between 1 µm and 10 µm are considered to be small particles, too small to be individually detected by a camera. Hence, the great amount of droplets within the pixel's solid angle change the intensity produced at a pixel. Volumetric scattering models such as attenuation and airlight [START_REF] Middleton | Vision through the Atmosphere[END_REF], [START_REF] Mccartney | Optics of the atmosphere: Scattering by molecules and particles[END_REF] are used to describe the effect of these small particles on the image.

In the last two decades, various methods [START_REF] Shi | Virtual scene simulation and implementation based on directx11[END_REF] have been developed for photo-realistic rendering of fog in computer graphics using scattering models and several studies [START_REF] Narasimhan | Vision and the atmosphere[END_REF], [START_REF] Garg | Vision and rain[END_REF] have been conducted to investigate the effect of bad weather on computer vision. Realistic fog simulation could allow evaluating the effects of foggy weather on computer vision algorithms under otherwise identical conditions (i.e. environment and camera setup).

In this document, we propose a method for the simulation of fog using realistic physical and photometric model for small particles and we investigate the influence of fog on computer vision algorithms. In the following we provide the details of our method for physical simulation of fog. We use the atmospheric scattering model to simulate fog and later show the advantage of heterogeneity for more realism. Finally, common features detectors are ran on images with simulated fog to evaluate the effects of fog.

II. SIMULATION OF FOG

In this section, we explain our method for rendering foggy images. As explained in [START_REF] Narasimhan | Vision and the atmosphere[END_REF], the small particles in the atmosphere creates a scattering medium for light. Light passing through a scattering medium is attenuated and distributed to other directions.

Authors are with the Robotics and Intelligent Transportation Systems (RITS) Team, INRIA Paris, 2 rue Simone Iff, 75012 Paris, France {raoul.de-charette, sule.kahraman} @inria.fr

A. Atmospheric Scattering Model

Atmospheric scattering model describes how suspended small particles alters the intensity and color of light due to attenuation and airlight.

First, we explain the effect of attenuation that fog creates. Due to attenuation, the radiance of a scene point decreases as its depth from the observer increases. Allards law [START_REF] Allard | Memoire sur lintensite et la portee des phares[END_REF] which is derived from the Bouguer's exponential law of attenuation describes this relation as follows:

E(d, λ ) = I o (λ )e -β (λ )d d 2 (1)
where d is the distance from the light source λ is the wavelength of light, I o is the radiant intensity of the point source, and β (λ ) is the total scattering coefficient of the medium.

Second effect of fog in the atmosphere is airlight. The latter causes the atmosphere to behave like a source of light due to the scattering of environmental illumination (e.i. sunlight, skylight, etc.) by particles in the atmosphere. Apparent brightness of a scene point increases with depth which means that airlight increases with path length, whereas attenuation causes scene radiance to decrease with path length. Total radiance at distance d from the camera, L(d, λ ) is

L(d, λ ) = L ∞ (1 -e -β (λ )d ) (2)
where L ∞ is the horizon radiance where the radiance of airlight is maximum when object at infinity. The total scattering coefficient is written as a function of λ in these equations because β depends on the ratio of droplet radius to wavelength and refractive index of water (which depends on temperature). However, for fog and dense haze, the influence of wavelength can be neglected [START_REF] Narasimhan | Vision and the atmosphere[END_REF]. This is because γ ≈ 0 in the Rayleighs law [START_REF] Middleton | Vision through the Atmosphere[END_REF], [START_REF] Mccartney | Optics of the atmosphere: Scattering by molecules and particles[END_REF], which describes the relationship between the scattering coefficient β , and the wavelength λ :

β (λ ) = Constant λ γ (3)
Hence, we assume the total scattering coefficient β (λ ) to be constant over because the spectral bandwidth of the camera is limited (visible light range for a gray-scale camera, and even narrower spectral bands when the camera is color) [START_REF] Narasimhan | Vision and the atmosphere[END_REF].

In order to physically simulate fog on images, we use the atmospheric scattering model which is the sum of the attenuation (eq. 1) and airlight models (eq. 2). Hence, the intensity of light received by the observer I(d, λ ) is the sum of directly transmitted light D(d, λ ) and airlight A(d, λ ) as described in [START_REF] Narasimhan | Vision and the atmosphere[END_REF]. Because the total scattering coefficient is independent of wavelength in the case of fog, we can describe the atmospheric scattering model as:

I(d) = D(d) + A(d) D(d) = I o e -β d A(d) = L ∞ (1 -e -β d ) (4)
When rendering fog on images, I(d) is the foggy image intensity as an RGB vector and D(d) is the product of transmission map t(d) and the clear scene image I o . Hence, I(d) is explicitly:

I(d) = I o e -β d + L ∞ (1 -e -β d ) (5)
Note that the spatially aligned depth map d is required. In virtual environment this is directly obtained from the rendering engine, whereas for real images we used MonoDepth [START_REF] Saxena | 3-d depth reconstruction from a single still image[END_REF] to predict depth from single RGB images.

B. Heterogeneous Fog

Using atmospheric scattering model in our renderings results smooth, evenly distributed homogeneous fog simulation, yet fog in the atmosphere is not homogeneous. Hence, images with homogeneous fog are not realistic. Common practice in the literature is to add 2D Perlin noise [START_REF] Perlin | Improving noise[END_REF] to create natural looking turbulence texture for rendering heterogeneous fog [START_REF] Zdrojewska | Real time rendering of heterogenous fog based on the graphics hardware acceleration[END_REF], [START_REF] Tarel | Vision enhancement in homogeneous and heterogeneous fog[END_REF]. However, 2D noise isn't realistic because fog itself is in a 3D space. 3D noise would better simulate spatial variations and potentially accurate spatio-temporal variations (assuming displacement is known). In our renderings we use Simplex noise simply because Simplex noise outperforms the more popular Perlin noise. Simplex noise, also developed by Ken Perlin, is designed to overcome the limitations of his classic Perlin noise algorithm. As S. Gustavson explains in [START_REF] Gustavson | Simplex noise demystified[END_REF], on the contrary to Perlin noise, Simplex noise has a gradient well-defined and continuous everywhere and it has a lower computational complexity. Although difficult to visualize, fig. 2 shows a binary representation of the simplex noise space.

We implemented the Simplex Noise [START_REF] Gustavson | Simplex noise demystified[END_REF] and substitute β in eq. 5, with β the noised scattering coefficient such that: β = N(β ). Where N(β ) is the noise function, randomly scaling β in the interval [β /2, 2β ].

The pinhole camera model [START_REF] Ikeuchi | Computer vision: A reference guide[END_REF] is used to retrieve the (X, Y, Z) camera coordinates from all image coordinates (x, y) and there associated depth (d). Camera coordinates are then mapped in the noise space to retrieve the ad-hoc noise value. Using this noise value to perturb the scattering coefficient would lead to acceptable simulation but non-realistic. Indeed it is often considered that the light attenuation is constant, but in fact all rays go through different volumes of atmosphere that exhibits varying scattering coefficients.

To make our model more precise, we average the scattering coefficients along the ray in each pixel. For each pixel, we calculate the noised scattering coefficient β for 10 discrete steps in [0, d] using the volume of 3D Simplex noise and average these scattering coefficients along the ray. This method is closer to the real life we are not assuming that the scattering coefficient is constant through the ray. In fig. 3, we render fog images with homogeneous (left) and heterogeneous (right) using the two aforementioned methods. As can be observed in fig. 3, fog with simplex noise is unevenly distributed and has different densities at different depths and points which produce more realistic images.

Figure 4 also illustrates the noise map of the heterogeneous fog shown above. Observing the noise maps, we see that the noise close to the camera is smooth and low frequency just like fog, yet in the sky area the frequency of fog is very high. This might be due to an error in our code or an unknown characteristic of the Simplex noise. Nevertheless, the high frequency noise in the sky should further be studied. In practice, note also that we map the camera space to the noise space with a scaling to ensure keeping a low frequency noise that resembles the density of fog. III. BENCHMARKING In this section we present and explain our results on the performance of some detection algorithms of Python's OpenCV library with the foggy conditions. After we rendered heterogeneous fog on images, we ran features detectors (Harris Corner, SIFT and SURF) algorithms with OpenCV and ran these algorithms on the foggy images. For Harris corner we compared the number of corners detected in the original clear scene image to the number of corners detected in foggy images with different scattering coefficient values ranging between 1 and 100. For SIFT and SURF, we compare the number of keypoints detected in the original clear scene image to the number of keypoints detected in the foggy images. As shown in Figure 5, performance of all these three detection algorithms decay very quickly with the influence of fog. IV. CONCLUSIONS With this work we have introduced a new method to simulate fog in images from RGB and depth data, accounting for physical and photometric models of small particles. Results are still very preliminary but acknowledge that fog has a great influence on the performance of basic features detectors which let foresee a probable influence on more complex computer vision algorithms. Future work include testing with object detection algorithms and investigating the presence of high frequency noise in distant areas (e.g. sky).

Fig. 1

 1 shows the original image along with a sample simulation of thick fog (β = 40). Additional visual results are shown in fig. 3. As the total scattering coefficient β increases, fog's density increases and the visibility decreases.

Fig. 1 .

 1 Fig. 1. Original image and sample homogeneous thick fog

Fig. 2 .

 2 Fig. 2. Binary visualization of the Simplex noise demonstration in 3D space

Fig. 3 .

 3 Fig. 3. Sample images showing homogeneous (left) vs. heterogeneous (right) fog, for different scattering coefficient.

Fig. 4 .

 4 Fig. 4. Noise map for the generation of heterogeneous fog

Fig. 5 .

 5 Fig. 5. Impact of fog on detection algorithms

ACKNOWLEDGMENT

This work has been conducted during an internship at Inria Paris under the supervision of Raoul de Charette.


Alexis Meyer

Raoul De Charette

Computing Ego Velocity from Scene Flow estimation

Keywords: Computer Vision, Image Flow, Ego motion, Intelligent Transportation System

This technical report describes a new approach that estimates the velocity of the ego vehicle using a monocular setup to compute an approximation of the scene flow.

I. INTRODUCTION

In the last two decades, several methods were proposed to estimate the ego motion of mobile robot as it is helps ego localization, trackings of objects and more generally improves scene understanding. Often ego motion extraction is done using LiDAR [START_REF] Zhang | Visual-lidar odometry and mapping: Low-drift, robust, and fast[END_REF], [START_REF] Zhang | Real-time depth enhanced monocular odometry[END_REF] sometimes with an optional IMU as in [START_REF] Scherer | River mapping from a flying robot: state estimation, river detection, and obstacle mapping[END_REF], or using vision [START_REF] Cvišić | Stereo odometry based on careful feature selection and tracking[END_REF], [START_REF] Buczko | Flow-decoupled normalized reprojection error for visual odometry[END_REF], [START_REF]How to distinguish inliers from outliers in visual odometry for high-speed automotive applications[END_REF]. For the latter, the fashion way is to use a stereo setup and estimate the displacement of the vehicle from the shift of keypoints in left/right images which provides a stereo disparity and thus RGB and Depth. As highlighted in [START_REF] Raudies | A review and evaluation of methods estimating ego-motion[END_REF] there is a large number of stereo techniques compared to the few monocular ones though some recent paper were proposed [START_REF] Mirabdollah | Fast techniques for monocular visual odometry[END_REF], [START_REF] Song | Parallel, real-time monocular visual odometry[END_REF], [START_REF] Song | Robust scale estimation in real-time monocular sfm for autonomous driving[END_REF].

We investigate in this document a method that requires only a monocular setup and compute the ego motion of the vehicle from the approximation of scene flow using a flat-world assumption. Our approach is novel and -though preliminary -the results obtained are promising. In the following we provide an overview of our methodology, details on our flow computation and the extraction of the optical flow. The document ends with our preliminary evaluations.

II. METHODOLOGY

We propose here a new approach to estimate ego velocity using a monocular camera setup. In contradiction to stereo, the task is made significantly harder since monocular vision can not provide any depth information.

To overcome the lack of depth data, we use a strong assumption: the ground is locally planar and all pixels below the horizon line are ground pixels. Using this assumption, we can compute the scene flow (velocity in the scene) from temporal image flow and calibration settings (section III). Although our assumption is strong and fails for many pixels (i.e. any non-ground pixels) we have found that the error induced by the latter is similar to a widely spread noise. Conversely, in the velocity/angle space the pixels validating our assumption (i.e. ground pixels) are well grouped and can be extracted through Gaussian fitting (section IV). Having done this the estimation of the ego motion is straightforward since the ground has a relative speed (i.e. in vehicle reference) equals to the inverse of the world vehicle velocity. We compare also Authors are with the Robotics and Intelligent Transportation Systems (RITS) Team, INRIA Paris, 2 rue Simone Iff, 75012 Paris, France {raoul.de-charette, alexis.meyer} @inria.fr different extraction methods: maximum, and various Mean Shift variations. In section V we evaluate our method on the widely used Kitti datasets [START_REF] Geiger | Are we ready for autonomous driving? the KITTI vision benchmark suite[END_REF] and discuss the promising results obtained.

III. TEMPORAL IMAGE FLOW COMPUTATION

The aim of this section is to explain how the temporal image flow has been exploited. The word temporal highlights the fact that this method only requires one single camera. The inputs of this method are the camera calibration parameters : internal parameters (focal length, magnification of the objective, optical center position) and the external ones (rotation and translation applied on the camera to pass from the vehicle to the camera coordinate system).

A. 2D image flow

Image flow computes the relative shift of pixel patches from image t to image t + 1, often through the minimization of a distance function in a restricting neighborhood.

In the literature there is two types of image flow, namely: sparse and dense. The former estimates shift of patches independently and fails for patches that could not be tracked (i.e. non-saliency patches). Dense image flow on the contrary interpolates patch motion to fill missing values. In our study, we used the sparse image flow Farneback method [START_REF] Farnebäck | Two-frame motion estimation based on polynomial expansion[END_REF] which approximates by a second degree polynomial, according to the following signal model:

f (x) ∼ x T Ax + b T x + c (1)
The identification of parameters A, b and c results of a weighted least squares approximation of the signal, where the weighting function normally is a Gaussian. The minimization of the image signal difference at time t and t + 1 in a given neighborhood provides an estimation of the temporal shift.

Figure 1 shows the temporal optical flow on the Kitti database.

The following parameters were used: winsize = 9, polygon n = 5, polygon sigma = 1.1, iteration number = 15 (refer to [START_REF] Farnebäck | Two-frame motion estimation based on polynomial expansion[END_REF] for parameters meaning). Note that to compensate the camera imaging noise a pre-smoothing is applied (9x9 Gaussian kernel).

In figure 1 the angle of the shifting was encoded as Hue whereas shift magnitude is encoded as Brightness. For example, bright red indicates a top right displacement whereas a dark . In the image flow, static elements are eccentric from the optical center (when the car is moving forward) whereas dynamic elements have an almost unpredictable behavior. Conversely, a fast overtaking car will induce a concentric image flow. As expected in such a dynamic scene the optical flow suffers from various problems, mainly: tracking precision due to a lack of saliency, tracking loss due border effects or large image displacement (more likely to occur with low camera frequency). In addition, the projection model implies smaller image shift for far-away displacement rather than close one, thus making harder to detect a dynamic object motion when far away.

The 2D image flow here obtained cannot be exploited as such. In fact, processing a 3D scene flow reveals to be the next step to a better understanding of the scene motion.

B. 3D scene flow

Forward projection geometry tells us that applying the projection matrix to a 3D vector coordinate point (X,Y, Z), we obtain its 2D imaged position (x, y):

  su sv s   =   f x s uv c u 0 0 f y c v 0 0 0 1 0       t x R 3x3 t y t z 0 0 0 1         X Y Z 1     (2)
and

  x y 1   =   su sv s   (3)
where u, v the image coordinates, s a scaling constant, f x , f y the focal length in pixels (respectively, width and height), c u , c v the coordinates of the optical center in the image, R 3x3 and t x , t y , t z the 3x3 rotation matrix and translation components to apply for camera to vehicle transformation. Using equation 2 we can compute a forward projection of a 3D point (scene to image), but we require backward projection (image to scene) which in practice can be done due to the loss of depth information in the projection process.

In spite of the unknown depth, we use inverse calibration projection geometry to compute the light ray direction for each pixel. The pseudo-inverse projection matrix M i is a matrix that if multipled from the right by the projection matrix is equal to the identity. If we multiply from left with M i , eq. 2 becomes:

M i   x y 1   =     X c Y c Z c 1     (4)
With eq. 4 we are able to obtain the projection vector (X c ,Y c , Z c ) of the pixel (x, y). Reader will note that because of the loss of s, there is no way to recover the depth along the vector. Instead our assumption is to assume that all pixels below the horizon belong to the ground surface which is approximated as locally planar. That is:

a * X + b * Y + c * Z = h (5)
With n = (a,b,c) the normal vector of the plane (in our case (0,1,0) because we consider that the flat ground is collinear to the plane Y=0). Hence, the pixel depth can be obtained from the intersection of the projection vector (X c ,Y c , Z c ) and the ground as expressed in eq. 5. This intersection occurs at (X,Y, Z), with:

     X = (X c -t x )d + t x Y = (Y c -t y )d + t y = h Z = (Z c -t z )d + t z (6)
where (X,Y, Z) are the real coordinates of a point of the scene, d the parameter of the system (the unknown variable) and h the altitude of the ground. Now that we are able to couple one pixel to a position in the scene we can turn the image displacement to scene displacement, in other words compute the scene flow from the image flow. Two components are needed to represent the velocity vector : the magnitude v of the velocity (in km/h) and its orientation angle θ .

     v = |d| * f * 3.6 |θ | = acos( V M .V n |V M ||V n | ) sign(θ ) = sign(V n X V M X -V n Z V M Z ) (7)
where d is the displacement vector of a considered element of the scene (of norm metric the meter) , V n the unit Z axis vector of the scene and f the frequency of the camera. Note that the multiplication 3.6 if for m/s to km/h conversion. Figure 2 shows an example output of our scene flow estimation, with the angle (θ ) encoded as Hue and velocity (v) encoded as Brightness. In the latter red (that is 0 • /360 • ) correspond to a forward relative motion in the scene, and cyan (that is 180 • ) is for backward relative scene motion. We can see that despite different optical flow our strong assumption leads to an accurate scene flow computation for all true ground pixels. This is visible in figure 2 where all road In the latter, the ground is represented with a unified color (dark cyan) which proves that our scene flow estimation is accurate for ground pixels.

pixels have similar dark cyan color (i.e. same scene velocity angle and magnitude). Interestingly reader will notice that even for non ground pixels that can not satisfy our assumption, the estimated scene flow is coherent. For example even if the black overtaking car is not the ground, the orientation angle of the velocity is coherent to the reality.

The 3D scene flow should now be analyzed in order to extract the ego motion of the vehicle. a) Mean Shift density estimator [START_REF] Fukunaga | The estimation of the gradient of a density function, with applications in pattern recognition[END_REF]. The principle of the latter is to converge iteratively towards the maximum of the distribution, assuming an underlying local Gaussian distribution. As for most density estimator the Mean Shift requires an initialization seed that is critical. In our case we used the maximum histogram as the initial seed. One of the main drawbacks of mean shift is the risk to extract a local maximum.

IV. EGO MOTION EXTRACTION

b) Mean Shift with multiple seeds was implemented to leverage the local maximum risk. Our 8 seeds our radially distributed on a circle centered on the maximum of the histogram and with radius equal to the uncertainty limit of the model. Our estimate here is the barycenter of the 8 meanshifts convergence.

c) Global Maximum was also used as a naive comparison.

We now discuss the evaluation of our complete methods and compare the different extraction approaches.

V. EVALUATION To evaluate the proposed method we used the Kitti benchmark database with the ad-hoc calibrations settings, and the associated velocity ground truth. Three datasets were evaluated from the training database 2 that are sequence 2, 8 and 9. Respectively refered as DB1, DB2, DB3. The complete extraction results of DB1 and DB3 are displayed in figure 4.

For all three databases we also detail the RMSE (Root Mean Square Error) which is a frequently used metric to measure the differences between two curves: here the ground truth and our extracted ego velocity.

Besides the complexity of the tested datasets (shadows, light saturation, lack of saliency, non-planar surfaces, traffic), the extracted velocities are coherent to the reality. When the car is turning, the estimated gaussian is curving to the opposite way as in fig. 6, which is logical since we observe the inverse ego motion vectors. Even when the car stops, the extracted ego motion reaches 0 km/h. Those comments underline the validity of our approach. The database DB 3 has 800 frames and many changes during the sequence (motion changes, lateral motion) whereas the two others are linear (no lateral motion, less changes) explaining its higher RMSE (DB 3 column in figure 5). In fact while the car meets changes, the image flow tracking algorithm performances can decrease even slightly causing the scene flow velocities to be wrongly estimated.

For high velocity, the pitch angle of the car changes abruptly (due to relief, pot holes, etc.) causing temporal inconsistencies locally in the scene flow extraction. As a consequence the extraction is more noisy for higher velocities.

In term of quantitative results, we were surprised to note that the maximum histogram reveals to be the best extraction among the three implemented. Nevertheless the mean shift multiple obtains close results, but the use of barycenter is probably not sufficient and better clustering should be applied for higher precision.

VI. CONCLUSIONS

This paper aimed to reveal an interesting and innovative way of ego motion extraction in the context of autonomous driving. The intentional use of a single camera setup made this approach challenging. Although the results are really preliminary and require additional work, they look promising and validate our strong assumption. Future works should consider using a more complex extraction method and a temporal filtering. Undoubtedly, this will improve performances.

  Fig. 1. Result of the sparse Farneback image flow [12] on Kitti dataset. Refer to text for details.

Fig. 2 .

 2 Fig.2. Image flow (top) and our corresponding estimation of the scene flow (bottom). In the latter, the ground is represented with a unified color (dark cyan) which proves that our scene flow estimation is accurate for ground pixels.

Fig. 3 .

 3 Fig. 3. Scene flow histogram : Discretisation step of 5 km/h and 5 degrees.

Figure 3

 3 Figure 3 shows the histogram of magnitude versus angle of the scene flow velocity. Since many pixels satisfy our assumption the histogram predominant distribution is located at the ground scene velocity location. We use histogram to extract the vehicle ego motion. To extract the velocity from angle/velocity distribution we compare three extraction methods: a) Mean Shift, b) Mean shift with multiple seeds, c) Global maximum.a) Mean Shift density estimator[START_REF] Fukunaga | The estimation of the gradient of a density function, with applications in pattern recognition[END_REF]. The principle of the latter is to converge iteratively towards the maximum of the distribution, assuming an underlying local Gaussian

Fig. 4 .

 4 Fig. 4. Ego motion extraction on DB1 and DB3 (cf. text for naming).

Fig. 5 .Fig. 6 .

 56 Fig. 5. RMSE calculation for each database. The maximum histogram seems to have the best results with in average 2.87 km/h.

ACKNOWLEDGMENT

This work has been conducted during an internship at INRIA Paris under the supervision of Raoul de Charette.