W Djema

L Giraldi

S Maslovskaya

O Bernard

Turnpike features in optimal selection of species represented by quota models

Keywords: Optimal control, turnpike properties, nonlinear systems, photobioreactor, strains selection, microalgae

The paper focuses on a generic optimal control problem (OCP) deriving from the competition between two microbial populations in continuous cultures. The competition for nutrients is reduced to a two-dimensional dynamical nonlinear-system that can be derived from classical quota models. We investigate an OCP that achieves species separation over a fixed time-window, suitable for a large class of empirical growth functions commonly used in quota models. Using Pontryagin's Maximum Principle (PMP), the optimal control strategy steering the model trajectories is fully characterized. Then, we provide sufficient conditions for the existence of a turnpike property associated with the optimal control and state-trajectories, as well as their respective costate trajectories. Indeed, we prove that for a sufficiently large time, the optimal strategy achieving strain separation remains most of the time exponentially close to an optimal steady-state defined from an associated simpler static-OCP. This turnpike feature is based on the hyperbolicity of the linearized Hamiltonian-system around the solution of the static-OCP. The obtained theoretical results are then illustrated on microalgae, described by the Droop model in dimension 5. The optimal strategy is numerically computed in Bocop (open source toolbox for optimal control ) with direct optimization methods.

Introduction

Natural selection drives the fate of the species that are permanently competing for limited survival sources. The individuals with a lower fitness in local conditions will eventually disappear from the ecosystem due to the competition for survival. Combined with mutations that may appear naturally over time, the Darwinian adaptation characterizing living organisms ensures that individuals which are better tailored to their local environment are continuously emerging and then progressively dominating their surroundings. In some instances, these evolution phenomena may be guided and redirected towards different fates. It is for instance possible to modify the outcome of species competition, leading to the emergence of enhanced organisms, by controlling and modulating the natural selection pressure [START_REF] Swenson | Artificial ecosystem selection[END_REF]). Such a selection process was then tamed for generating naturally more efficient microogranisms [START_REF] Almeida | Selection and application of microorganisms to improve oil recovery[END_REF]). Within this scope, some microalgae strains Corresponding author: Walid Djema. have been recently improved by setting-up a Darwinian selection pressure within artificial cultures such as photobioreactors (see, e.g., [START_REF] Gachelin | Enhancing PUFA-rich polar lipids in Tisochrysis lutea using adaptive laboratory evolution (ALE) with oscillating thermal stress[END_REF]; [START_REF] Huertas | Estimat-ing the capability of different phytoplankton groups to adapt to contamination: herbicides will affect phytoplankton species differently[END_REF]; [START_REF] Bonnefond | Continuous selection pressure to improve temperature acclimation of Tisochrysis lutea[END_REF]).

Microalgae attract considerable attention since they have a great potential to produce less impacting biofuels [START_REF] Chisti | Biodiesel from microalgae[END_REF]; [START_REF] Wijffels | An outlook on microalgal biofuels[END_REF]). In addition to reducing atmospheric CO 2 , they are also considered as promising elemental bricks for green chemistry and profitable new protein sources. Moreover, they are even used to improve wastewater treatment (see also, [START_REF] Walsh | Algal food and fuel coproduction can mitigate greenhouse gas emissions while improving land and water-use efficiency[END_REF]; [START_REF] Odjadjare | Potential biotechnological application of microalgae: a critical review[END_REF]). However, despite their large interest, microalgae are still in their infancy and hitherto mainly grown for cosmetic or pharmaceutical industries. Their exploitation on a larger scale, particularly for biofuel production, is conditioned by the overcoming of various technical challenges, including the necessity of enhancing the strain productivity since mainly wild species are currently cultivated. This explains the worldwide race during the last decade for selecting more suitable microalgae strains (see, e.g., [START_REF] Aravantinou | Selection of microalgae for wastewater treatment and potential lipids production[END_REF]; [START_REF] Rodolfi | Microalgae for oil: Strain selection, induction of lipid synthesis and outdoor mass cultivation in a low-cost photobioreactor[END_REF]; [START_REF] Yoo | Selection of microalgae for lipid production under high levels carbon dioxide[END_REF]).

In this paper, the objective is to develop an effective control strategy ensuring that the most suitable strains will dominate the continuous artificial culture in finite time. For instance, in biofuel production, desirable microalgae strains can be those which grow faster, with higher lipid content, and which are also well adapted to local environment and seasons.

We recall that when n distinct strains coexist in a chemostat with limited substrate, a well-known principle -the competitive exclusion principle CEP [START_REF] Smith | The theory of the chemostat: dynamics of microbial competition[END_REF])-predicts that the species requiring the smallest amount of nutrients to grow at a rate equal to the dilution rate wins the competition, while the other species will be asymptotically washed out. However, depending on the culture conditions, the winner of the competition is not in general necessarily the one with the highest lipid content or unrivalled productivity. Therefore, our goal is to determine an optimal control strategy so that the species of interest is the one that dominates in finite time. To achieve our objectives, we use an approach that is based on optimal control theory (Pontryagin's principle, [START_REF] Pontryagin | Mathematical theory of optimal processes[END_REF]) and turnpike features (see, e.g., [START_REF] Trélat | The turnpike property in finitedimensional nonlinear optimal control[END_REF]). There are other works for optimal selection of microorganisms in minimal-time [START_REF] Bayen | Optimization of the separation of two species in a chemostat[END_REF], 2017)), based on the simpler Monod-type dynamics describing the growth of bacteria. The microalgae behavior is however more complex since the growth of these microorganisms does not depend on the substrate concentration in the photobioreactor (like in the Monod model), but rather on the quantity of nutrients stored inside the cells [START_REF] Bernard | Modelling of microalgae culture systems with applications to control and optimization[END_REF]; [START_REF] Bernard | Hurdles and challenges for modelling and control of microalgae for CO 2 mitigation and biofuel production[END_REF]). The Droop model [START_REF] Droop | Vitamin B 12 and marine ecology. IV. The kinetics of uptake growth and inhibition in Monochrysis lutheri[END_REF]) is widely used to represent algal growth, since it takes into account the additional dynamics governing the internal cell-quota storage [START_REF] Smith | The theory of the chemostat: dynamics of microbial competition[END_REF]; [START_REF] Droop | 25 years of algal growth kinetics, a personal view[END_REF]). Here we use a general form of quota models, extending the classical Droop's formulation to include different forms of growth rate. In particular, we suppose in this work that an ad'hoc controller is able to track the desired nutrient concentration in the medium. As a consequence, the classical 5 dimensional quota model (studied for instance in Djema et al. (2020a,b)) reduces to a system of dimension 2. In addition, this work generalizes the approach in [START_REF] Djema | An Optimal Control Strategy Separating Two Species of Microalgae in Photobioreactors[END_REF](Djema et al. ( , 2020b)))), not only because it considers a generic class of systems relevant for microalgae (using different growth rates), but also because there are no assumption on the initial conditions, which are often assumed to be already confined to an attractive positively invariant set [START_REF] Bayen | Optimization of the separation of two species in a chemostat[END_REF]; [START_REF] Djema | An Optimal Control Strategy Separating Two Species of Microalgae in Photobioreactors[END_REF]Djema et al. ( , 2020b))). The main novelty of our work is to prove the existence of a turnpike property characterizing the optimal solutions. The turnpike has been recently revisited in the optimal control literature and it is attracting more and more attention (see, e.g., [START_REF] Trélat | Integral and measure-turnpike properties for infinite-dimensional optimal control systems[END_REF]; [START_REF] Porretta | Long time versus steady state optimal control[END_REF]; [START_REF] Zaslavski | Turnpike theory of continuous-time linear optimal control problems[END_REF]; [START_REF] Trélat | Integral and measure-turnpike properties for infinite-dimensional optimal control systems[END_REF]; [START_REF] Lohéac | Minimal controllability time for finite-dimensional control systems under state constraints[END_REF]; [START_REF] Gugat | On the turnpike phenomenon for optimal boundary control problems with hyperbolic systems[END_REF]). The turnpike property states that the optimal solutions consist of three pieces: two transient short-time arcs at the boundaries of the time interval, and a long-time behavior remaining exponentially close to a static solution [START_REF] Trélat | The turnpike property in finitedimensional nonlinear optimal control[END_REF]). From a practical standpoint, turnpike features are also advantageous for deriving sub-optimal controls that are easy to implement on real industrial systems.

The paper is organized as follows. Sect. 2 is devoted to the statement of the general optimal control problem (OCP), formulated for a generic two-dimensional dynamical system, which includes competition between algal species. In Sect. 3 we state the main theorem dealing with the general OCP. Its proof is provided in Sect. 4, which relies on Pontryagin's maximum principle (PMP, [START_REF] Pontryagin | Mathematical theory of optimal processes[END_REF]) and also on the existence of a turnpike property along the optimal solutions. In Sect. 5, we focus on an application involving a variable yield model that includes the classical Droop's model. Thus, based on the results in Sect. 3-4, we solve the problem of microalgae selection under several growth-rates. Finally, in Sect. 5.5, the optimal control is determined in different cases and the turnpike feature is illustrated using a direct optimization method implemented in Bocop (optimal control toolbox, [START_REF] Bonnans | BOCOP: an open source toolbox for optimal control -A collection of examples[END_REF]).

Optimal control problem (OCP)

We consider a generic nonlinear two-dimensional dynamical system, suitable for competition between two populations in which the dynamics are coupled through the feeding control-flow. Typically, the competition between two species (e.g., phytoplankton, bacteria, yeasts,...) for a limiting substrate (w), whose concentration is controlled in a continuous bioreactor, can be represented as follows:

ẏ1 = -f 1 (y 1 )y 1 + g 1 (w), ẏ2 = -f 2 (y 2 )y 2 + g 2 (w), (1)
where the f i (•) are regular nonnegative real-valued and increasing functions representing the growth rate of species i with respect to the intracellular limiting nutrient y i . The functions g i (•) are also regular nonnegative real-valued and increasing standing for the flux of limiting nutrient into the cell. System (1) typically appears as a subsystem of quota models, where y i is the internal quota of the i-th species, i.e. the intracellular amount of the element for which the involved species are competing. The Droop's model is extensively studied later in Sect. 5 as the most famous example of this class of systems. However, the model (1) can also be found in more complicated metabolic dynamical systems [START_REF] Wiechert | Modeling and simulation: tools for metabolic engineering[END_REF]) where the intracellular element y i is diluted by growth. The resulting system (1) is positive, i.e., for positive initial conditions the trajectories remain positive.

In general, the dynamics of the system (1) is at higher dimension (5 in the case of Droop's model), but here we assume that the limiting substrate has been controlled by an ad'hoc controller tracking the concentration w(t). Indeed, the way this subsystem can be extracted from the full model is developed in Sect. 5 in the typical case of Droop's model describing microalgae. Without loss of generality, we consider throughout this paper that species 1 is the one of interest and must become predominant in the population in finite time. Our objective is to investigate the control strategies for this competition between two initially co-existing populations by controlling the nonnegative and bounded input w over a fixed time-interval [0, T f ].

Statement of the general OCP

The set of admissible controls associated with system (1) is defined by w ∈ W, s.t., W is a subset of L ∞ loc (R + ), the space of locally integrable functions on every compact set on R + . More precisely, we set,

w ∈ W = {w : [0, T f ] → [0, w m ] | w(•) ∈ L ∞ loc (R + )}
, where w m is a sufficiently large positive constant (the maximum influent substrate concentration). Then, we define now the main OCP that reads,

       Maximize w∈W C [0,T f ] = T f 0 f 1 (y 1 ( )) -f 2 (y 2 ( ))d ,
where, y = (y 1 , y 2 ), y(•) is the solution of (1) corresponding to the control w(•).

(

In what follows, the terminal conditions of the states y i , i = 1, 2, are assumed to be free and their initial conditions are positive. Besides, let us define (y * i , w * ), i = 1, 2, the steady state of (1), where, y * i > 0 and w * ∈ [0, w m ]. It follows that, f i (y * i )y * i = g i (w * ), for i = 1, 2, and since the map, y i → f i (y i )y i , is increasing, we get, w * = δ i (y * ), where, δ i (y) = g -1 i (f i (y i )y i ). Thus, y * i (w * ) = δ -1 i (w * ), i = 1, 2, for all fixed w * ∈ [0, w m ]. Then now, we define for all admissible y * i > 0 and w * ∈ [0, w m ] the function,

∆(w * ) = f 1 (δ -1 1 (w * )) -f 2 (δ -1 2 (w * )). ( 3
)
Throughout this work, we assume that the maximizer of ∆ is the unique point w ∈ (0, w m ), as illustrated in Fig. 1. The cases where w is zero or w m are pointless and can be trivially addressed outside this framework. Let us remark that under these considered assumptions, OCP

(2) has a maximizer using Filippov Theorem (see, e.g., [START_REF] Cesari | Optimization -Theory and applications[END_REF]). In fact, the generic form of the functions f i (δ -1 i (•)), involving one intersection point over (0, w m ) as in Fig. 1, ensures the existence of control functions maximizing and minimizing the criterion C [0,T f ] , as a direct consequence of the competitive exclusion principle (CEP), (see, [START_REF] Smith | The theory of the chemostat: dynamics of microbial competition[END_REF]). i (w)), in their typical general form studied throughout this work, where both situations

f1(δ -1 1 (w)) > f2(δ -1 2 (w)) and f1(δ -1 1 (w)) < f2(δ -1 2 (w)) hold over [0, wm].
The resulting function ∆ in (d) satisfies the hypothesis w ∈ (0, wm).

Pontryagin's maximum principle (PMP)

The Hamiltonian associated with OCP ( 2) is defined by,

H =Φ(λ y1 , λ y2 , w) -λ y1 f 1 (y 1 )y 1 -λ y2 f 2 (y 2 )y 2 + λ 0 [f 1 (y 1 ) -f 2 (y 2 )] , (4) where,
Φ(λ y1 , λ y2 , w) = λ y1 g 1 (w) + λ y2 g 2 (w).

(5)

It is classical to consider that λ 0 = 1 since OCP ( 2) is stated as a maximization problem. The PMP implies the existence of absolutely-continuous λ yi , which are the costate of y i , satisfying, λyi = -∂H/∂y i , for i = 1, 2. In addition, since the final states y i (T f ) are free, it follows that the transversality conditions are given by,

λ yi (T f ) = 0, for i = 1, 2. (6)
Finally, the PMP implies that the control function

w ∈ W satisfies, for almost all t ∈ [0, T f ], max w(t)∈[0,wm] H(y i (t), λ yi (t), w(t)). (7)

Introducing the static-OCP

The static-OCP (see, e.g., [START_REF] Trélat | The turnpike property in finitedimensional nonlinear optimal control[END_REF]) associated with the OCP (2), for the system at steady state is:

   Maximize w∈[0,wm] ∆(w) = f 1 (y 1 (w)) -f 2 (y 2 (w)), under the constraints f i (y i )y i = g i (w), i = 1, 2. (8)
In other words, the static-OCP seeks to maximize the function ∆ defined in (3) when system (1) is at steady state. Under the assumption that w ∈ (0, w m ) is the unique point maximizing ∆, as in Fig. 1, the definition of the static-OCP is extended to include the co-state variables using the Hamiltonian equation. Indeed, according to the Lagrange multipliers rule, there exists, (λ y1 , λ y2 ) ∈ R 2 , s.t., the solution (y 1 , y 2 , w), implies that y 1 , y 2 , λ y1 , λ y2 , w also satisfies:

       ∂H ∂λy i (y 1 , y 2 , λ y1 , λ y2 , -1, w) = 0, -∂H ∂yi (y 1 , y 2 , λ y1 , λ y2 , -1, w) = 0, ∂H ∂w (y 1 , y 2 , λ y1 , λ y2 , -1, w) = 0. ( 9
)
More precisely, using H given in (4) we deduce that the steady-state satisfying the static-OCP (8) is solution of,

             -f i (y i )y i + g i (w) = 0, i = 1, 2, λ y1 (f 1 (y 1 ) + f 1 (y 1 )y 1 ) -f 1 (y 1 ) = 0, λ y2 (f 2 (y 2 ) + f 2 (y 2 )y 2 ) + f 2 (y 2 ) = 0, g 1 (w)λ y1 + g 2 (w)λ y2 = 0, (10)
In the general case illustrated in Fig. 1, we define the unique static optimal steady-state that satisfies (10), denoted by,

E static = w, y 1 , y 2 , λ y1 , λ y2 . (11)
In Sect. 5, an explicit determination of E static is given for typical f i and g i functions in quota models.

3 Turnpike main result : Theorem 1

Let us consider the perturbed trajectories around the static optimal steady-state E static defined in (11). More precisely, we perform the classical change of coordinates,

δy i (t) = y i (t) -y i , δλ yi (t) = λ yi (t) -λ yi , for i = 1, 2,
and we consider the shifted control, δw(t) = w(t) -w. Thus, as in [START_REF] Trélat | The turnpike property in finitedimensional nonlinear optimal control[END_REF]), we can check that, δw(t) = -H -1 ww (H wy δy(t) + H wλ δλ y (t)), where, δy = (δy 1 , δy 2 ) T and δλ y = (δλ y1 , δλ y2 ) T . In what follows, we set z = (δy, δλ y ) T . Using the Hamiltonian (4), we notice that

H yw = ∂ 2 H ∂y∂w E static = (0, 0),
where,

H ww = ∂ 2 H ∂w 2 E static = λ y1 g 1 (w) + λ y2 g 2 (w), H yy = diag{ -λ y1 (2f 1 (y 1 ) + y 1 f 1 (y 1 )) + f 1 (y 1 ), -λ y2 (2f 2 (y 2 ) + y 2 f 2 (y 2 )) -f 2 (y 2 )}, H λy = diag {-f 1 (y 1 ) -y 1 f 1 (y 1 ), -f 2 (y 2 ) -y 2 f 2 (y 2 )}, H wλ = (g 1 (w), g 2 (w))
, and thus we end up with,

H λw H -1 ww H wλ =   g 2 1 (w) λy 1 g 1 (w)+λy 2 g 2 (w) g 1 (w)g 2 (w) λy 1 g 1 (w)+λy 2 g 2 (w) g 1 (w)g 2 (w) λy 1 g 1 (w)+λy 2 g 2 (w) g 2 2 (w) λy 1 g 1 (w)+λy 2 g 2 (w)   .
Let us define the following matrices and constants,

-H λw H -1 ww H wλ H yy = b 1 b 2 b 3 b 4 , b 23 = b 2 b 3 , a i = -∂(fi(yi)yi) ∂yi E static , H λy = diag{a 1 , a 2 }.
Consequently, we can show that the eigenvalues κ j , j = 1, . . . , 4, of the matrix Ψ defined in ( 14) are given by,

κ 1 = 1 √ 2 a 2 1 + a 2 2 -b 1 -b 4 + ς 1/2 , κ 2 = -κ 1 , κ 3 = 1 √ 2 a 2 1 + a 2 2 -b 1 -b 4 -ς 1/2 , κ 4 = -κ 3 , (15)
where,

ς = (a 2 1 -a 2 2 ) -(b 1 -b 2 ) 2 -4b 23 .
Now, we are ready to state the main result of the paper.

Theorem 1 Let T f > 0 be a large fixed final-time and E static be the unique global solution of the static-OCP (8), s.t. w ∈ (0, w m ) maximizes ∆ in (3).

x The solution of OCP (2) with fixed initial conditions

y i (0) ∈ R * + , i = 1, 2, correspond to the singular control function w s ∈ W, s.t., w s (t) = {0, w m } , for all t ∈ [0, T f ], (16)
i.e., the optimal control maximizing H in (4) is neither a bang 0 nor a bang w m .

y Moreover, if a 1 = a 2 , H ww < 0, and,

Re(κ j ) = 0, for j = 1, . . . , 4, (17)
where κ j are the eigenvalues of the matrix Ψ given in (15), then there exist positive constants , c 1 , c 2 and a time-instant

T * > 0, s.t., for all T f > T * , if, 2 i=1 y i (0) -y i + λ yi < , (18)
then the optimal solutions of OCP (2)

satisfy ∀t ∈ [0, T f ], 2 i=1 y i (t) -y i + λ yi (t) -λ yi + w s (t) -w ≤ c 1 e -c2t + e -c2(T f -t) , (19)
where c 1 does not depend on T f and the optimal control is the singular w s in x.

Remark 1 Assuming that T f is sufficiently large is classical in turnpike framework [START_REF] Trélat | The turnpike property in finitedimensional nonlinear optimal control[END_REF]).

The proof of Theorem 1 provides an estimate on the lower-bound of T f . This bound is satisfied in biological applications, e.g., metabolite production (Yegorov et [START_REF] Trélat | The turnpike property in finitedimensional nonlinear optimal control[END_REF]. Notice that the main theorem from [START_REF] Trélat | The turnpike property in finitedimensional nonlinear optimal control[END_REF] cannot be applied to our class of optimal control problems. Firstly, the control is bounded in the considered class of optimal control problem whereas it is unbounded in [START_REF] Trélat | The turnpike property in finitedimensional nonlinear optimal control[END_REF].

Secondly, the matrix -H yy + H yw H -1 ww H wy is not positive definite in our case. For instance, choosing f i (y i ) =

aiyi-bi yi

, it follows that 2f 1 (y 1 ) + y 1 f 1 (y 1 ) = 0, and then we deduce that, -H yy + H yw H -1 ww H wy = -H yy = diag{f 1 (y 1 ), -f 2 (y 2 )}, which does not satisfy the positive definite condition in [START_REF] Trélat | The turnpike property in finitedimensional nonlinear optimal control[END_REF]. This point is illustrated numerically in Sect. 5.5 using Droop's functions. Therefore, the hypotheses of the second part of Theorem 1 are slightly different from the hypotheses of the main theorem in [START_REF] Trélat | The turnpike property in finitedimensional nonlinear optimal control[END_REF].

The proof of Theorem 1 is divided in two parts. The first step concerns the statement x and derives from the PMP. The second step is based on the hyperbolicity of the linearized system ż = Ψz + Λ(z), where z = (δy, δλy) T . The arguments similar to those in [START_REF] Pighin | Long time behaviour of Optimal Control problems and the Turnpike Property[END_REF]; [START_REF] Trélat | The turnpike property in finitedimensional nonlinear optimal control[END_REF] are sketched, the proof is detailed in the associated technical report: [START_REF] Djema | Turnpike features in optimal selection of species represented by quota models: extended proofs[END_REF].

4 Proof of Theorem 1

x Using (4)-( 5), we note that the PMP implies that the admissible control function satisfies, for almost all t ∈ [0, T f ], the maximization condition, max

w(t)∈[0,wm] H(y i , λ yi , w) = max w(t)∈[0,wm] Φ(λ yi , w). ( 20
)
At first sight, we deduce from (20) that any extremal solution should be a concatenation of bang and singular arcs over [0, T f ]. Indeed, using (5) we notice that in the case where λ y1 (t) ≥ 0 and λ y2 (t) ≥ 0, the control satisfying (20) is the bang w(t) = w m . Similarly, if λ y1 (t) ≤ 0 and λ y2 (t) ≤ 0, then it follows that the control satisfying (20) is the bang w(t) = 0. And finally, when λ y1 (t)λ y2 (t) < 0 at any t ∈ [0, T f ], we deduce from the fact that g i is a positive increasing function that the control satisfying ( 20) is the singular control denoted by w s (t), which maximizes Φ and which is neither 0 nor w m . However, using the co-states dynamics given by,

λy1 = λ y1 (f 1 (y 1 ) + f 1 (y 1 )y 1 ) -f 1 (y 1 ), λy2 = λ y2 (f 2 (y 2 ) + f 2 (y 2 )y 2 ) + f 2 (y 2 ), (21)
where f i are continuously increasing, we deduce that,

λ y1 = f 1 (y 1 ) f1(y 1 )+f 1 (y 1 )y 1 > 0 and λ y2 = -f 2 (y 2 ) f2(y 2 )+f 2 (y 2 )y 2 < 0.
Then, since λ yi is an equilibrium point of ( 21) and using the transversality conditions (6), we deduce that,

0 < λ y1 (t) < λ y1 , λ y2 < λ y2 (t) < 0, (22)
for all t ∈ [0, T f ), where the latter inequalities derive from the fact that ∂(f i (y i )y i )/∂y i > 0 and f i (y i ) > 0 in (21). In particular, (22) imply that for all t ∈ [0, T f ),

λ y1 (t)λ y2 (t) < 0. ( 23
)
Consequently, we conclude from (23) that the control w s (t) is the one that satisfies the PMP and thus steers over the time window [0, T f ] the trajectories from y 1 (0), y 2 (0) to y 1 (T f ), y 2 (T f ). This concludes the proof of x.

y Now, let us focus on the eigenvalues of the matrix Ψ defined in (14). Firstly, let us observe that,

det(Ψ -κI 4 ) = det H λy -κI 2 -H λw H -1 ww H wλ -H yy -H yλ -κI 2
, where I p is the p × p identity matrix. We notice that H yy and -H yλ -κI 2 are diagonal and thus they are commutative. Therefore, it follows that,

det(Ψ -κI 4 ) = det (H λy -κI 2 )(-H λy -κI 2 ) -H λw H -1 ww H wλ H yy .
Consequently, det(Ψ -κI 4 ) = 0, leads to,

(a 2 1 -κ 2 -b 1 )(a 2 2 -κ 2 -b 4 ) -b 23 = 0, ( 24
)
where the constants a's and b's are defined before Theorem 1. From the equation ( 24) we readily get the eigenvalues κ j , j = 1, . . . , 4, given in (15). Under the suitable assumptions expressed in (17), we ensure that the studied optimality system is hyperbolic, i.e., all eigenvalues of Ψ have nonzero real parts. Next, by Theorem 6 of [START_REF] Molinari | The time-invariant linear-quadratic optimal control problem[END_REF], it follows that there exist symmetric 2 × 2 matrices P -, P + maximal and minimal solutions of the Riccati equation, H yλ P + P H λy -H yy -P H λw H -1 ww H wλ P = 0. Moreover, the corresponding matrices Ω + = H λy -H λw H -1 ww H wλ P + and Ω -= H λy -H λw H -1 ww H wλ P - are such that eigenvalues of Ω + have strictly positive real part and eigenvalues of Ω -have strictly negative real part, i.e., Ω + is unstable and Ω -is stable.

We denote by, V := I I

P -P +
, where I is the 2 × 2 identity matrix. Straightforward calculations show that Ψ = V -1 ΨV = diag{Ω -, Ω + }. We denote by z the new variable s.t. V z = z, and then it follows from (12) that,

   ż1 (t) = Ω -z1 (t) + θ 1 (z 1 (t), z2 (t)) , ż2 (t) = Ω + z2 (t) + θ 2 (z 1 (t), z2 (t)) , (25)
where the θ i , i = 1, 2, satisfy the o-little property as in (13). At this juncture, we point out that under different hypothesis, we have rewritten the optimality conditions in the form (25) similarly to [START_REF] Trélat | The turnpike property in finitedimensional nonlinear optimal control[END_REF] and [START_REF] Pighin | Long time behaviour of Optimal Control problems and the Turnpike Property[END_REF]. It can be shown (see Lemmas 2-3, Djema et al. (2021)), that there exists r > 0 s.t. for all (z 0 1 , zf 2 ) satisfying, z0 1 = z1 (0), zf 1 = z2 (T f ), (z 0 1 , zf 2 ) ∈ B(0, r), the corresponding solution (z 1 (t, (z 0 1 , zf 2 ), z2 (t, (z 0 1 , zf 2 )) of system (25) admits for any

t ∈ [0, T f ] the bound,    z1 (t) ≤ C κ z0 1 e -κt + e -κ(T f -t) zf 2 Θ 1 ( z1 c 0 ) z2 (t) ≤ C κ zf 2 e -κ(T f -t) + e -κt z0 1 Θ 2 ( z2 c 0 ) ( 26
)
Next step is to verify that there exist T * > 0 and > 0 in (18) s.t., for any T f > T * and y(0) satisfying ( 18), the initial and final conditions (z 0 1 , zf 2 ) satisfy,

z0 1 + z2 (0) = y 0 -y, P -z1 (T f ) + P + zf 2 = -λ y . (27)
For that, we use the fixed point argument as in [START_REF] Pighin | Long time behaviour of Optimal Control problems and the Turnpike Property[END_REF], by defining the map F : B(0, r) → R 4 ,s.t.,

F z0 1 zf 2 = y 0 -y -z2 (0) -P -1 + (P -z1 (T f ) + λ y )
, and showing that F admits the unique fixed point (z 0 1 , zf 2 ). This step is detailed in Appendix D.4 in [START_REF] Djema | Turnpike features in optimal selection of species represented by quota models: extended proofs[END_REF]. It is worth mentioning that at this stage of the proof, we also get the lower-bound T * of T f and from (18). Finally, using (26) for T f > T * , we obtain for all t ∈ [0, T f ],

δy(t) = z1 (t, (z 0 1 , zf 2 )) + z2 (t, (z 0 1 , zf 2 )) ≤ C κ z0 1 e -κt + e -κ(T f -t) zf 2 Θ 1 (r) +C κ zf 2 e -κ(T f -t) + e -κt z0 1 Θ 2 (r) .
Up to defining a bigger constant C κ , we deduce that,

δy(t) ≤ C κ e -κt + e -κ(T f -t) , ∀t ∈ [0, T f ]. (28)
Note that C κ is independent of T f . Similarly, we get

δλ y (t) ≤ C κ e -κt + e -κ(T f -t) , ∀t ∈ [0, T f ]. ( 29
)
Note that there exists R = R(r) s.t. δy + δλ y ≤ R.

Then, using additional calculations for δw(t), we get,

δw(t) = w(y, λ y ) -w ≤ sup B((ȳ, λy),R) ∇ y w δy + sup B((ȳ, λy),R) ∇ λy w δλ y , (30)
for all t ∈ [0, T f ]. Combining ( 28), ( 29) and ( 30) leads us to the exponential turnpike estimate given in ( 19) and thus concludes the proof of Theorem 1.

5 Application to a higher dimensional system

Algal competition using Droop's model

In this section, we focus on a model of competition between two microalgae populations in dimension 5, with the aim of achieving species separation over a fixed time window. We tweak this issue using a system reduction in order to rewrite the optimization objective in term of the generic OCP (2). Let x 1 and x 2 be the biomass concentrations of two different strains of microalgae, initially coexisting in a bioreactor with one limiting substrate (s), s.t., s(t) ∈ [0, s in ], for all t ≥ 0, where s in is the constant substrate influent concentration. Each microalgae species uptakes the free nutrient s. The intracellular amount of this limiting element is the cell quota q i , with i = 1, 2.

A general class of quota models, including the classical Droop's model [START_REF] Droop | Vitamin B 12 and marine ecology. IV. The kinetics of uptake growth and inhibition in Monochrysis lutheri[END_REF]), writes,

       ṡ = (s in -s)D - 2 i=1 ρ i (s)x i , qi = ρ i (s) -µ i (q i )q i , ẋi = [µ i (q i ) -D] x i , (31)
The function µ i in (31) stands for the growth rate of the i-th species, while ρ i is its nutrient uptake rate. The dilution rate D and the substrate input s in can generally be manipulated. The functions ρ i and µ i are nonnegative and increasing bounded functions, s.t., 0 ≤ ρ i (s) ≤ ρ mi , and, 0 ≤ µ i (q i ) ≤ µ mi , where ρ mi and µ mi are positive constants. Typically in Droop's model, the uptake rate ρ i (s) is expressed in terms of Michaelis-Menten kinetics,

ρ i (s) = ρ mi s K si + s , ( 32
)
where K si is a positive constant. The most classical form -known as Droop's growth-for µ i is given for i = 1, 2 by,

µ i (q i ) = µ i∞ 1 - Q i min q i , q i ≥ Q i min . ( 33
)
The minimum intracellular quota Q i min > 0 defines the quota under which cell division cannot occur. Other forms have been used to describe the growth rate such as the Caperon-Meyer equation [START_REF] Caperon | Nitrogen-limited growth of marine phytoplankton -II. Uptake kinetics and their role in nutrient limited growth of phytoplankton[END_REF]),

µ i (q i ) = µ i∞ (q i -Q i min ) q i -Q i min + k qi , q i ≥ Q i min , k iq > 0, ( 34
)
another possible form that we also use in this work is,

µ i (q i ) = µ i∞ q ni i q ni i + k ni qi , n i > 1. ( 35
)
Notice that the choices (34)-( 35) lead to nonlinear terms µ i (q i )q i in (31). In all these cases, the resulting quota is bounded for all t ≥ 0, i.e., Q i min ≤ q i (t) ≤ q mi , where q mi is the maximum internal storage rate [START_REF] Bernard | Hurdles and challenges for modelling and control of microalgae for CO 2 mitigation and biofuel production[END_REF]). In the sequel, the quota variables q i will stand for y i in model ( 1) and we will prove that: (i) OCP (2) achieves microalgae separation, through maximization of relative concentration x 1 (T f )/x 2 (T f ), (ii) Theorem 1 applies for different functions µ i (given by ( 33)-( 35)), i.e., turnpike-type behaviors characterize microalgae selection under various growth rates.

Controlling the substrate concentration

Now, let us assume that the free substrate concentration s is controlled by an ad'hoc controller.For instance, by playing with the concentration of the influent substrate s in and with the dilution rate D, the concentration s can follow a desired (smooth enough) reference trajectory s(t) = s * (t). Different control designs can be used such as linearizing control (Bastin (1990)), asymptotic control [START_REF] Mailleret | Nonlinear adaptive control for bioreactors with unknown kinetics[END_REF]) or other approaches that demonstrated their efficiency on microalgae (see, e.g., [START_REF] Bernard | Modelling of microalgae culture systems with applications to control and optimization[END_REF]; [START_REF] Tebbani | CO2 Biofixation by Microalgae: Modeling, estimation and control[END_REF]). By assuming a perfect tracking of the substrate towards w(t), the system (31) then reduces to,

qi = ρ i (w) -µ i (q i )q i , ẋi = [µ i (q i ) -D] x i , (36)
This system is triangular, and the dynamics of x i results from the dynamics of q i . Consequently, we can focus exclusively on the quota-dynamics q i , for i = 1, 2.

Maximizing the relative biomass density

Species separation is achieved when maximizing the relative density of the species of interest x 1 at T f , i.e., finding the optimal control maximizing the quantity,

C(T f ) = x 1 (T f ) x 2 (T f ) . ( 37
)
Let us state the following key result to end up with a functional that does not depend on the states x i , but still maximizes the relative concentration at the final-time.

Proposition 1 Any solution of the class of systems (36) maximising C [0,T f ] defined in OCP (2) is equivalently maximizing CT f in (37), where the functions f i in OCP

(2) stand for the growth rates µ i in system (36).

Proof: Indeed, the cost function in OCP (2) writes,

C [0,T f ] = T f 0 [µ 1 (q 1 (t)) -D] -[µ 2 (q 2 (t)) -D] dt. Using the state dynamics, µ i (q i ) -D = ẋi xi , we get, C [0,T f ] = [ln(x 1 (t)) -ln(x 2 (t))] T f 0 = ln x1(T f )x2(0) x1(0)x2(T f ) . Thus, maximizing C [0,T f ] indeed consists in maximizing CT f . The advantage of considering the criterion C [0,T f ]
instead of the natural criterion CT f is that it is independent of the dynamics of the states x i , for i = 1, 2. Thus, the optimization problem can now be reduced to an OCP as in (2), involving the cost C [0,T f ] and associated with the reduced-system (in the form (1)),

   q1 = -µ 1 (q 1 )q 1 + ρ 1 (w), q2 = -µ 2 (q 2 )q 2 + ρ 2 (w), (38)
where the growth rate µ i is given by one of the functions (33)-( 35), the uptake rate ρ i is given by (32), and the control w is equivalent to the substrate s in (31) once regulated by the ad'hoc controller. This is especially true since x in Theorem 1 states that the solution of the OCP (2) corresponds to the singular control function with no bangs over the time-window [0, T f ], and consequently the optimal control w # (t) has thus the required smoothness to be tracked in the original system (31) through the ad'hoc controller. To sum up, we can see that finding a control for system (31) approaching as close as possible the solution w # (t) of the OCP studied in Sections 2-3 will provide an excellent suboptimal control for the OCP of the ful-system (31) associated with the criterion (37).

Optimal synthesis of species separation

Now, we are ready to apply the results obtained in Sect. 2-3 to the problem of species selection. Firstly, the static-OCP associated with (38) and

C [0,T f ] states,    Maximize w∈[0,wm] ∆(w) = µ 1 (q 1 (w)) -µ 2 (q 2 (w)),
where -µ i (q i )q i + ρ i (w) = 0 for i = 1, 2 , (39) under the steady-state constraints -µ i (q i )q i +ρ i (w) = 0 for i = 1, 2. Similarly to Section 2, the solution of (39) -including the co-states-is denoted by, E static = q 1 , q 2 , λ q1 , λ q2 , w . This optimal steady-state satisfies ( 9)-( 10), using this time the quota-coordinates in (38), under the standard assumption of Section 2, i.e., w ∈ (0, w m ) as in Fig. 1, which is generally satisfied for ρ i and µ i in ( 32)-( 35), as illustrated later in numerical simulations. A case in point is the classical Droop's model associated with ( 32)-( 33), in which case the function ∆ rewrites,

∆(w) = µ 1 (δ -1 1 (w)) -µ 2 (δ -1 2 (w)), (40)
where, δ i ( ) = ρ -1 i (µ i ( ) ), and since, ρ

-1 i ( ) = Ki ρim-, for all ∈ [0, ρ im ), we obtain, δ i ( ) = K i -Qi min σi-
, where, σ i = ρim µi∞ + Q i min . This leads to, δ -1 i ( ) = σi +KiQi min +Ki , and consequently we end up in the Droop's case with,

µ i (δ -1 i ( )) = ρ im σ i + K i Q i min , i = 1, 2. ( 41
)
Using (41) we can compute successive derivatives of ∆(w) in the Droop's case, i.e., ∆ (w) = ρ1mK1Q1 min (K1Q1 min +σ1w) 2 -ρ2mK2Q2 min (K2Q2 min +σ2w) 2 , and, ∆ (w) = κ 1 (w)κ 2 (w), where, κ i (w) = -2Kiµ 2 i∞ Qi min ρim(µi∞Qi min +ρim) (Kiµi∞Qi min +(µi∞Qi min+ρim)w) 3 . Thus, by looking for the points that satisfy ∆ (w) = 0 and ∆ (w) < 0 in the case where µ i is given by (33), we deduce that the steady-state E static = w, q 1 , q 2 , λ q1 , λ q2 , satisfies,

q i = 1 µi∞ (ρ i (w) + µ i∞ Q i min ) , for i = 1, 2, λ q1 = -Q1 min q 2 1
, and, λ q2 = Q2 min q 2 2 . Numerically, we check that only one numerical value w belongs to (0, w m ) for a sufficiently large w m , leading then to the configuration in Fig. 1 (Fig. 4 in Droop's case). We refer to Appendices A-B in [START_REF] Djema | Turnpike features in optimal selection of species represented by quota models: extended proofs[END_REF] for more details on the practical determination of E static in the Droop's case. In a similar way, we can determine the analytic expression of E static for different growth-rates ((34)-( 35)). Therefore, under suitable conditions a unique E static exists and satisfies w ∈ (0, w m ), and thus in light of Theorem 1 we derive the following statement as a direct consequence.

Claim 1 (Application of Theorem 1) Let fix any biological parameters for the functions ρ i in (32) and µ i in (33)-( 35), i = 1, 2, in system (31), s.t., there exists a unique E static with w ∈ (0, w m ), solution of the static-OCP (39) then, x the optimal control w # (t) for all t ∈ [0, T f ] is singular, i.e., it does not involve bang arcs 0 and w m over [0, T f ]. y Moreover, if the corresponding conditions (17) and (18) are fulfilled, the singular optimal control (substrate concentration) w # (t) and its related quota q i (t), as well as their respective co-states λ qi (t), i = 1, 2, exhibit over the fixed time-horizon [0, T f ] an exponential turnpike behavior characterized in (19).

Note that the overlapping of favorable features (namely the facts that w solution of the static-OCP is remote from the bounds 0 and w m , the optimal control is inherently singular, i.e., w s (t) ∈ (0, w m ) for all t ∈ [0, T f ], and also the exponential turnpike keeping w s (t) almost constant most of the time) results in strong smoothness of the optimal control and the optimal trajectories, as illustrated in the next section. Nonetheless, this smoothness does not imply that the resulting ad'hoc controller of system (31) should be smooth as well, as shown in Djema et al. (2020b) where the obtained dilution-based control is of type bang-singular. To sum up the methodology provided in this paper: the first step is to determine the optimal point E static s.t. w ∈ (0, w m ), solution of the static-OCP corresponding to a given dataset of biological parameters and functions in the system (31). Next, we compute the matrix Ψ defined around E static as in ( 14). Then, the conditions stated in y-Theorem 1 ensure the existence of a turnpike-type behavior of the optimal solution of the reduced system (1). The numerical computation of the eigenvalues κ j , j = 1, . . . , 4, which are defined in (15), is illustrated in different cases of the functions µ i in Section 5.5. Finally, it is worth noting that the conditions (17) in Theorem 1 appear to be systematically satisfied in numerical simulations, leading accordingly to a turnpike behavior in all the evaluated cases.

Numerical synthesis through direct methods

The direct-optimization approach performed in this section illustrates the results of Theorem 1 applied to the quota-model in Sect. 5 (System (38)) for the different growth functions ((33)-( 35)) illustrated in Fig. 2. In this perspective, the state variables and the control function of system (38) are discretized, by setting a subdivision of the time-horizon [0, T f ] with fixed T f , and the substrate-control w is discretized over W. The state dynamics q 1 and q 2 (38) is discretized with Lobatto scheme (Betts (2010)), i.e. RK method of type Lobatto-IIIC, based on an implicit trapezoidal rule, of order 6. suitable for stiff-nonlinear problems. Thus, numerical direct methods transform OCP (2) (stated in Sect. 2, associated with System (38)) through the discretization step, into a nonlinear programming problem in finitedimension (see [START_REF] Djema | Turnpike features in optimal selection of species represented by quota models: extended proofs[END_REF], Sect. 5.5). Using the numerical values of the biological parameters in Fig. 2, Fig. 2. Various functions used in Droop and Caperon-Meyer models [START_REF] Caperon | Nitrogen-limited growth of marine phytoplankton -II. Uptake kinetics and their role in nutrient limited growth of phytoplankton[END_REF]; [START_REF] Droop | Vitamin B 12 and marine ecology. IV. The kinetics of uptake growth and inhibition in Monochrysis lutheri[END_REF]). From left to right: absorption functions ρi in (32) with, ρ1m = 9.1 (µmol.L -1 .d -1 ), K1 = 0.08 (µmol.L -1 ) ρ2m = 9.6 (µmol.L -1 .d -1 ), K2 = 0.25 (µmol.L -1 ). Growth rate µi in (33) with, µ1∞ = 1.6 day -1 , Q1 min = 0.95 µmol/L, µ2∞ = 2.2 day -1 , Q2 min = 1.5µmol/L. Growth rate µi in (34) with, µ1∞ = 0.75day -1 , Q1 min = 1.6µmol/L, kq1 = 0.9µmol/L, µ2∞ = 0.9day -1 , Q2 min = 2µmol/L, kq2 = 0.6µmol/L. Growth rate µi in (35) with, µ1∞ = 1.5day -1 , n1 = 4, kq1 = 3µmol/L, µ2∞ = 1.15day we performed several direct optimizations in Bocop with appropriate settings (tolerance for NLP solver: 10 -34 , 2000 time-steps). No issue is encountered regarding the initialization and starting points for solving the studied problem. In all the following examples, the conditions of Theorem 1 are satisfied, and the optimal solutions satisfy the turnpike property (19) as stated in Claim 1.

Example 1 We set T f = 15 days and we consider that the growth functions µ i are in the Droop's form (33). In this case, the sufficient conditions given in [START_REF] Trélat | The turnpike property in finitedimensional nonlinear optimal control[END_REF] are not satisfied (see Remark 3). The function ∆

(w) = µ 1 (δ -1 1 (w)) -µ 2 (δ -1 2 (w)) is given in Fig. ( 3
). In addition, the static steady-state E static and the eigenvalues of the matrix Ψ are given in Tab. . Notice also that in this case H yy = diag{0.2997, -1.0443}. According to Theorem 1 (see Claim 1), the turnpike property (19) holds as shown in Fig. 4.

Example 2 We set T f = 30 days and we consider that the growth functions µ i are in the form (34), with the biological parameters given in Fig. 2. Results are in Tab. 1. The turnpike holds according to Theorem 1, as described in Claim 1 and illustrated in Fig. 5. Notice that the sufficient conditions for the turnpike existence given in [START_REF] Trélat | The turnpike property in finitedimensional nonlinear optimal control[END_REF] are not satisfied, since H yy = diag{0.116993, -0.939631}.

Example 3 We set T f = 50 days. The growth functions µ i are in the form (35), with the biologi-cal parameters given in Fig. 2. Results are in Tab. 1. Similarly to the previous examples, here we get, H yy = diag{0.183898, -0.210028}. However, the turnpike characteristics hold in the optimal solution according to Theorem 1 and Claim 1, as shown in Fig. 6.

In all these examples (Fig. 456), the optimal solutions of the OCP (2), settled in relatively large time [0, T f ], consist of three pieces: the first and the last ones are relatively short arcs, and the middle piece being a relatively long-time arc staying close to the optimal steady-state solution of the associated static-OCP (8).

Conclusion

The optimal control has been fully-characterized through a complete study combining Pontryagin's principle and turnpike frameworks relying upon classical hyperbolic system properties. The turnpike theory [START_REF] Trélat | The turnpike property in finitedimensional nonlinear optimal control[END_REF]; [START_REF] Pighin | Long time behaviour of Optimal Control problems and the Turnpike Property[END_REF]) considerably simplifies the understanding of the control strategy and provides a clear insight into the optimal control structure. By investigating the OCP and its associated static version, the main result in this paper proves that the optimal control is always singular and provides sufficient conditions ensuring the existence of turnpike behaviors characterizing the optimal trajectories. We had to upgrade the existing turnpike results so that they can apply to the specificity of our models. The considered class of systems naturally appears in the framework of competition between two microbial or algal species, but a similar structure is likely to play a role for other classes of problems. In general, solving an OCP in dimension 5 is a challenging issue (see Djema et al. (2020a,b)). Thus, the microalgae optimal selection problem was reduced in this paper considering the core dynamics in dimension 2. This reduction approach using a back-stepping philosophy might also prove to be efficient in other contexts (see, e.g., [START_REF] Yegorov | Optimal control of bacterial growth for the maximization of metabolite production[END_REF]; [START_REF] Caillau | The turnpike property in maximization of microbial metabolite production[END_REF], and this will be the main focus of our future work.

  Fig. 1. (a-b) provide the generic form of the functions fi and gi for i = 1, 2 (with arbitrary units), studied in this paper. (c) illustrates the resulting fi(δ -1i (w)), in their typical general form studied throughout this work, where both situations f1(δ -1 1 (w)) > f2(δ -1 2 (w)) and f1(δ -1 1 (w)) < f2(δ -1 2 (w)) hold over[0, wm]. The resulting function ∆ in (d) satisfies the hypothesis w ∈ (0, wm).

Fig. 3 .

 3 Fig.3. The functions µi(δ -1 i (w)) and ∆(w) using the data given in Fig.2with µi as in (33).

Fig. 4 .Fig. 5 .Fig. 6 .

 456 Fig.4. The optimal control and optimal trajectories given by Bocop in Example 1.

  and thus around the unique point E static we end up with,

							and,
							Ψ =	H λy -H λw H -1 ww H wλ	,
							-H yy	-H yλ
	δ ẏ(t) δ λy (t)	= Ψ	δλ y (t) δy(t)	+	Λ 2 (δy(t), δλ y (t)) Λ 1 (δy(t), δλ y (t))	, (12)
	for all t ∈ [0, T f ], where Λ i satisfies,
		Λ (δy,δλy) → 0	0.	(13)

i (δy, δλ y ) (δy, δλ y ) --------→

  Values of the static steady-state E static and eigenvalues of the matrix Ψ defined at E static for the 3 examples.

				-1 ,
	n2 = 3, kq2 = 2 µmol/L.
		Example 1	Example 2	Example 3
	w	0.021	0.015	0.074
	q 1	2.164	0.015	3.91
	q 2	1.84	0.015	2.69
	λq 1	0.44	0.322	0.173
	λq 2	-0.20	-0.1377	2.690
	κ1,2	2.332 ±0.81i	1.74 ±1.09i	±3.02
	κ3,4 -2.332 ±0.81i -1.74 ±1.09i	±0.82
	Table 1

Acknowledgements: This work was supported by the ADEME BioMSA and Maximic ANR projects in France. The authors


Vincenzo Caglioti

Pierluigi Taddei

Planar Motion Estimation Using an Uncalibrated General Camera

Generic camera odometry can be performed using two images of the same plane, such as the ground plane even in the case of non central cameras. It is possible to recover both the angle and rotation center describing a generic planar motion on the ground plane if the center is visible in both images. We present an algorithm to recover these two parameters from an initial set of correspondences and, furthermore, to estimate the motion flow related to any point on the ground plane. In this situation the motion flows are given by a set of "concentric" closed curves around the rotation center. By considering two subsequent ground plane motions and their related motion flows, we show that it is possible to perform a rectification of the plane up to a scale factor. We provide experimental results which validate our approach.

Introduction

Robot localization is a fundamental process in mobile robotics application. One way to determine the displacements and measure the movement of a mobile robot is dead reckoning systems. However these systems are not reliable since they provide noisy measurements and tend to diverge after few steps [START_REF] Borenstein | Measurement and correction of systematic odometry errors in mobile robots[END_REF]. Visual odometry, i.e. methods based on visual estimation of the motion through images captured by one or more cameras, is exploited to obtain more reliable estimates. Furthermore, using an omni-directional camera increases the performance of the visual data, thus reducing the error cumulation effect.

Usually, 3D reconstruction from images taken by a moving uncalibrated perspective camera go through auto-calibration. Auto-calibration from planar scenes requires either non-planar motion [START_REF] Triggs | Autocalibration from planar scenes[END_REF], or several planar motions with different attitudes of the camera w.r.t. the ground plane [START_REF] Knight | Linear auto-calibration for ground plane motion[END_REF]. In the case of standard uncalibrated perspective camera, instead, it as been shown that visual odometry is feasible from a single planar motion [START_REF] Caglioti | Uncalibrated visual odometry for ground plane motion without auto-calibration[END_REF].

Generic cameras, which associate one projection ray to each individual pixel, represent the most general mathematical model to describe a projection system [START_REF] Grossberg | A general imaging model and a method for finding its parameters[END_REF]. In this case the relation among points on the ground plane and the image plane is not trivial and thus standard visual odometry techniques cannot be applied. Moreover, due to the difficulty of the calibration step, requiring a calibrated general camera is a strong assumption which often could not be assumed as possible.

In this paper we perform visual odometry using general cameras by exploiting a single planar motion on a textured ground floor. Our algorithm can be applied to an arbitrary projection model without requiring the camera to be central. Notice that the initial and final position induced by any planar motion are always related by a particular roto-translation. Apart pure translations, the motion can described by a center of rotation and a rotation angle. We estimate the motion flows induced by the roto-translation, i.e. the set of image points that matches an initial point at some stage of the motion [START_REF] Ramalingam | Towards generic self-calibration of central cameras[END_REF]. For such motions, in fact, the flow curves are represented by closed paths around the imaged center of rotation. Since the camera is general, though, the curves are not given by ellipses, as in the perspective case, but have a more general form. A sampling of the motion flows is recovered from the initial set of correspondences either by fusion, when more than two images of the same motion are available, or by interpolation. We follow this latter approach since we assume that only two images are available for the same roto-translation on the ground plane.

Given the sampled motion flows, we show how to estimate with good precision both the imaged center and angle of rotation. The center represents the only invariant point of the function describing the motion flows, i.e. its path coincide with a point. The rotation angle, instead, is recovered by counting the number of samples required to cover a closed path around the center. We then show that, given the previous estimations, it is possible to recover a continuous function representing the motion flows. This function describes the path of all image points for any rotation angle.

Finally, given multiple subsequent planar motions, we show that it is possible to perform a rectification of the ground plane up to a scale factor. Any closed path, in fact, represents the image of a circumference centered on one of the rotation centers. We consider different flows related to the various motions and their intersections on the image plane. These points represent the intersections of the respective circumferences on the ground plane. It is possible to obtain the metric reconstruction by constraining the angles swept by the motion flows from one intersection to the other on each closed path.

In the following sections elements on the ground plane are represented by capital letters whereas elements on the image plane are represented by lower case letters. We use the dot notation Γ (p, •) to indicate that the function Γ (p, α) is considered as if the parameter α was fixed, i.e. Γ (p, •) depends only on p.

The paper is organized as follows. In Section 3 we formulate the problem and present the center and rotation angle estimation algorithms in Section 4. Section 5 shows how to recover a consistent motion flow using a trivariate b-spline. Section 6 describes how to restore the ground plane exploiting two different planar motions, i.e. given two set of motion flows around different image points. In Section 7 we show good results on real data related to the visual odometry and preliminary results on the ground plane rectification. Finally, Section 8 concludes the paper.

Previous Work

In the last years methods to estimate the robot motion (ego-motion) based on visual information provided by cameras have gained attention. Early methods were based on estimation of the optical flow from image sequence in order to retrieve ego-motion. McCarthy and Barnes [START_REF] Mccarthy | Performance of optical flow techniques for indoor navigation with a mobile robot[END_REF] presented a review of the most promising methods.

Other approaches exploited stereo vision. Nister et al. [START_REF] Nister | Visual odometry[END_REF] proposed a method based on triangulation between stereo pairs and feature tracking in time-sequence of stereo pairs, without any prior knowledge or assumption about the motion and the environment. Takaoka et al. [START_REF] Takaoka | 3d map building for a humanoid robot by using visual odometry[END_REF] developed a visual odometry system for a humanoid robot based on feature tracking and depth estimate using stereo pairs. Agrawal and Konolige [START_REF] Agrawal | Real-time localization in outdoor environments using stereo vision and inexpensive gps[END_REF] proposed an integrated, real-time system involving both stereo estimate in the disparity space and a GPS sensor in a Kalman filter framework. GPS-based systems can be sufficiently accurate for large areas but it can be not used in indoor environments and require a support framework, which prevent their use, e.g., for planetary exploration. For such application, Mars Exploration Rover [START_REF] Cheng | Visual odometry on the mars exploration rovers -a tool to ensure accurate driving and science imaging[END_REF] employed a features detection in a stereo image pair that are tracked from one frame to the next; using maximum likelihood estimation the change in position and attitude for two or more pairs of stereo images is determined.

Davison [START_REF] Davison | Real-time simultaneous localization and mapping with a single camera[END_REF] proposed a real-time framework for ego-motion estimation for a single camera moving through general unknown environments. The method was based on a Bayesian framework that detect and track a set of features (usually corner, points or lines). Assuming the rigidity in the scene, the feature image motion allows to estimate the motion of the camera; therefore the complete camera trajectory and a 3D map of all the observed features can be recovered. Wang et al. [START_REF] Wang | Visual odometry based on locally planar ground assumption[END_REF] and Benhimane and Malis [START_REF] Benhimane | Homography-based 2D visual servoing[END_REF] estimated the homographies between image sequences to retrieve robot motion but both methods require camera calibration. Caglioti and Gasparini [START_REF] Caglioti | Uncalibrated visual odometry for ground plane motion without auto-calibration[END_REF], instead, performed visual odometry without going through calibration.

Visual odometry system based on non perspective cameras have also been proposed. Gluckman and Nayar extendend standard techniques to omnidirectional cameras [START_REF] Gluckman | Ego-motion and omnidirectional cameras[END_REF]. Bunschoten and Krose [START_REF] Bunschoten | Visual odometry from an omnidirectional vision system[END_REF] used a central catadioptric camera to estimated the relative pose relationship from corresponding points in two panoramic images via the epipolar geometry; the scale of the movement is subsequently estimated via the homography relating planar perspective images of the ground plane. Corke et al. [START_REF] Corke | Omnidirectional visual odometry for a planetary rover[END_REF] developed a visual odometry system for planetary rover based on a catadioptric camera; they proposed a method based on robust optical flow estimate from salient visual features tracked between pairs of images. Ramalingam et al. [START_REF] Ramalingam | Towards generic self-calibration of central cameras[END_REF][START_REF] Ramalingam | Minimal solutions for generic imaging models[END_REF] used (motion) flow lines to calibrate a general central camera from general 3D motions whereas our work deals with planar motions but is applicable also to non-central cameras.

Problem Statement

Given an initial image I and a final image I taken by a camera which undergoes a planar motion w.r.t. the visible ground plane, we aim at recovering the angle ᾱ and the imaged center of rotation c of the camera. We will consider as given a set of correspondences {p → p } between the two images related to points on the ground plane.

A motion flow is represented by a function relating any image point p to an image curve containing all points that match p at some stage during the motion. In the case of a planar motion w.r.t. the ground plane, the motion flow of any imaged point is repre-sented by a closed curve around the center of rotation [START_REF] Ramalingam | Towards generic self-calibration of central cameras[END_REF]. These curves are represented by conics in the case of perspective cameras, but have more complex form in the case of generic projection (see Figure 1). In order to estimate ᾱ and c we will recover a function

Γ ᾱ : R 2 × Z → R 2
which represents a sampling of the continuous motion flows induced by a planar rototranslation around c.

Given an image point p and a integer value n, Γ ᾱ(p, n) represents the matching point after n rotations of magnitude ᾱ around c.

Notice that Γ ᾱ(p, 0) represents an identity while Γ ᾱ(p, 1) relates points in image I with points in image I and, conversely, Γ ᾱ(p, -1) relates points in image I with points in image I.

Let p 0 and p 1 refer to two matching image points, i.e. observing the same 3D point in I and I . Let p 2 refer to the point that in I matches to point p 1 in I. Similarly let p 3 be the point that in I matches to point p 2 in I, and so forth (see Figure 1). This sequence may be expressed using function Γ ᾱ, and in particular:

p 0 = Γ ᾱ(p 0 , 0) p 1 = Γ ᾱ(p 0 , 1) p 2 = Γ ᾱ(p 0 , 2) = Γ ᾱ(p 1 , 1) = Γ ᾱ(Γ ᾱ(p 0 , 1), 1) . . . p n = Γ ᾱ(p 0 , n) = Γ ᾱ(. . . Γ ᾱ(p 0 , 1), . . . , 1)
n times which means that Γ ᾱ(p, n) may be evaluated by applying n times Γ ᾱ(•, 1) to point p.

In order to estimate Γ ᾱ we thus estimate Γ ᾱ(•, 1) (and by symmetry Γ ᾱ(•, -1)) using the known set of correspondences. We show that this approach is valid in order to recover ᾱ and c although the solution suffers from drifting, i.e. the function is likely to be wrong for large values of n. In fact the recovered sequence of matching points Γ ᾱ(p, •) does not represent the sampling of a closed curve. In Section 5 we show how to handle this problem in order to recover motion flows while imposing that they are closed curves. This step is essential in order to perform a rectification of the ground plane.

Planar Motion Estimation

In order to estimate Γ ᾱ(•, 1) we exploit the set of correspondences {p → p } between I and I . The two images are registered using a bivariate b-spline function Γ . This mathematical model allows one to express the registration problem as a sparse linear least square minimization.

In particular we aim at minimizing the reprojection error of feature points in I plus an additional smoothing term, which depends on the second derivatives of the b-spline.

Γ ᾱ = arg min (E r [Γ ] + λ s E s [Γ ]) (1)
where

E r [Γ ] = i |Γ (p i , 1) -p i | 2 (2)
E s [Γ ] = j |Γ uu (q j , 1)| 2 + |Γ vv (q j , 1)| 2 + 2 |Γ uv (q j , 1)| 2 . (3)
Here Γ uu , Γ vv and Γ uv represent the second partial derivatives of Γ , which are sampled at points {q j }, and λ s is a weighting constant.

In order to deal with catadioptirc cameras and in general with cameras with large field of view we exploit a b-spline in polar coordinates in order to better fit the image region which should represent the ground plane.

Imaged Center Estimation

The only invariant point of the function

Γ ᾱ(•, 1) is point c, in fact ∀k ∈ Z Γ (c, k) = c. ( 4
)
and all other points are displaced according to the motion. Their flows are represented by closed curves around c. In particular these sequences of points all rotate around c either clockwise or counter-clockwise Consider a line crossing c. In the general case, all points on one half line are displaced in one half-space and all points on the other half line are displaced in the opposite half space due to the previous property (see Figure 2). This observation leads to simple, although effective, way to compute the imaged center.

We start by considering an initial segment crossing a point c 0 . The segment is sampled at n sites and the displacement given by Γ ᾱ(•, 1) is evaluated. We then detect a point c 1 which is the point on the segment that divides the sites in two group w.r.t. the displacement orientation. We then perform the same step on an orthogonal segment passing through c 1 detecting point c 2 and so forth.

By iterating the process, while reducing the length of the segment used, the sequence of {c i } converges to c.

Angle Estimation

We estimate ᾱ using the following procedure. Let us consider a generic image point p and the sequence of points given by Γ ᾱ(p, .).

In particular let us consider a reference line b crossing both c and the chosen initial point p (see Figure 3).

Suppose that the recovered function Γ ᾱ is correct and that the rotation angle is an integral fraction f of 2π. In this case point Γ ᾱ(p, f ) coincides with p, i.e. the sequence has a period of f . Each application of Γ ᾱ(•, 1) represents a rotation of 2π f . In the case that the rotation angle is not an integral fraction of 2π, it is estimated by considering the least common multiple of 2π and α. This is practically evaluated by considering a large sequence of points which performs several loops around c. By counting the number of times the sequence traverses the line b and dividing by the number of points used we obtain a good estimate of ᾱ.

Notice however that our initial estimate of Γ ᾱ(•, 1) is not accurate and suffers from drifting (see Figure 3). We observed experimentally that this simple method accurately estimates the rotation angle.

Motion Flow Estimation

In the previous section we showed that, generally, the recovered function Γ ᾱ suffers from drifting. The most notable effect of this problem is that the sequence of points representing a motion flow does not lie on a closed curve. Moreover this function is a sampling of the real flows induced by the planar motion.

To perform the ground plane rectification we need a function which is constituted by closed curves and, moreover, which is continuous w.r.t. the angle α. For this reason we will exploit Γ ᾱ(•, 1) to recover a continuous function

H : R 3 → R 2 (5)
such that H(p, α) is the matching point of p after a rotation of α around the center C.

In this formulation H has the following properties

H(p, 0) = p ∀p ∈ R 2 , (6) H(c, α) = c ∀α ∈ R, (7)
H(p, α) = H(p, α + 2π) ∀p ∈ R 2 , ∀α ∈ R, (8)
H(•, ᾱ) = Γ ᾱ(•, 1), (9) H(•, -ᾱ) = Γ ᾱ(•, -1). ( 10
)
We estimate function H using a trilinear periodic b-spline. In particular we perform an optimizaton enforcing the previous properties, that can all be expressed as linear constraints w.r.t. the b-spline parameters, and an additional smoothing term. The optimization is implemented as a sparse linear least square minimization.

Ground Plane Rectification

Once H is estimated, we can associate to any image point a continuous closed curve. This curve represents the image of a circumference on the ground floor. In particular all the closed curves given by H are the images of a set of concentric circumferences around the center of rotation C.

Consider a camera performing a first planar motion where the initial and final positions are given by images I 1 and I 2 . Consider a second planar motion with a different center of rotation, where the inital and final position are given by images I 3 and I 4 . Using the technique described in Section 5 we obtain two different continuous functions H 1 and H 2 . These functions continuously describe two rotations around the two centers.

Let C 1 and C 2 be the centers of rotation on the ground floor of the two motions and c 1 and c 2 their imaged points on one of the images. Moreover let L be the line passing through C 1 and C 2 . Given any point P on the ground plane there exist two circumferences centered in C 1 and C 2 intersecting in P . These circumferences intersect also in a symmetric point Q w.r.t. to line L. Point Q is the matching point of P after a rotation of angle α 1 , where α 1 is the angle formed by P and Q w.r.t. the first circumference. Similarly α 2 is the angle related to the other circumference (see Figure 4). Consider a point p on the image plane, and its corresponding point P on the ground plane. The point q corresponding to the symmetric point Q of P lies on the same closed curves, H 1 (p, •) and H 2 (p, •), passing thorough p since these two curves must represent the images of the two ground plane circumferences. By detecting the intersections of these two curves it is possible to recover both point q and the angles α 1 and α 2 .

Without loss of generality, consider a reference system for the ground plane such that C 1 lies at the origin and C 2 at point (1, 0). In this case line L corresponds to the x axis. Consider a pair of points P and Q symmetric w.r.t. L. Line P C 1 forms with L an angle α1 2 whereas line QC 1 forms an angle -α1 2 . Similarly line P C 2 forms an angle 2π-α2 2

and line QC 2 an angle -2π-α2

2

. Moreover the points

H 1 (p, α1 2 ), H 1 (p, π + α1 2 ), H 2 (p, α2
2 ) and H 2 (p, π + α2 2 ) represent the imaged intersections of the two circumferences with L.

Thus for any point p on the image plane we can collect 8 constraints, 4 constraining points to lie on line L and 4 constraining the position of P and Q. These constraints are used to perform a rectification of the ground plane by sampling the image plane using N points and recovering the 8N constraints. The rectification is solved registering a bivariate b-spline using a sparse linear least square optimization involving a smoother terms and the collected constraints.

Experimental Validation

We tested our algorithm on real images taken with a central camera, although the algorithm is more general and could be applied to non central cameras as well. We manually extracted a set of correspondences and then recovered function Γ ᾱ(•, 1). The ground truth for the angle and center of rotation was roughly measured during the experiments setup. We run the algorithm to recover the imaged center c starting from different initial configurations. In all test cases the recovered solution was the same, showing its independence from the initial starting solution. We, then, let the number of parameters of the bivariate b-spline vary, i.e. increasing the allowed amount of deformation. In particular we noticed that the recovered solution tends to be closer to the measured ground truth as the number of parameters increases but as this number increases further the estimation tends to diverge by several pixels. This is probably due to overfitting.

In our experiments with 70 correspondences the estimation was good up to a bspline using a grid of 35 × 35 knots. In all experiments the standard deviation of the recovered points was less than 6 pixel on 1024 × 768 images (Figure 5 shows the recovered centers on two of such experiments).

The iterative algorithm for the angle recovery converged rapidly to the average value in all our experiments. An example of the evaluation is shown in Figure 6. The initial point used for the angle estimation slightly influences the recovered value. We ran the algorithm several times from initial random solutions. In all experiments, the standard deviation was less than 0.2 degree (see Figure 7).

As for the center of rotation, we ran the algorithm several times varying the number of parameters of the b-spline. The estimation was robust since the standard deviation was always below 0.4 degree. Also, for this experiment, the accuracy increased with the Above this threshold the estimation diverged (see Figure 7).

Besides our experiments on the motion estimation (rotation center and angle), we also performed preliminary experiments on the ground plane rectification. Given two different rotations, we evaluated the functions H 1 and H 2 (such as the ones shown in Figure 8).

We applied the algorithm for the ground floor rectification using these two functions and obtaining the result shown in Figure 9. The regions containing more extracted features were correctly reconstructed. The rectification, though, presented some errors where the features are too sparse, in particular at the center of the image, due to the occlusion of the camera, and at the outer regions, due to the resolution of the sensor. We are currently working on improvements of the ground plane rectification through the use of multiple planar motions.

Conclusion

We presented a method to perform visual odometry using general cameras for planar motions when the equivalent roto-translation center is visible. In particular given the initial and final image of a planar motion on a ground plane we showed how to recover both the imaged center and the angle of rotation.

Moreover we presented an algorithm to estimate a continuous function describing the motion flows of any point on the image plane. By using two of these function, i.e. by considering two subsequent planar motions, we showed how to perform a rectification of the imaged ground plane.

The angle and center estimation showed to be robust to the number of degrees of freedom used in the modelization, i.e. the parameters of the b-spline function, and to the initial solution.

The mathematical model used to describe H, although composed by closed curves, does not constrain those to be concentric, i.e. not to intersect each other. Moreover there are few constraints to describe the behavior of this function at different values of α. A possible solution could be fusing point correspondences of multiple images related to the same planar motion. Another possibility is exploiting more than two planar motions to increase the number of constraints for the rectification. Future works will be aimed at improving the estimation of such a function. The final rectification, in fact, depends on how well the function H describes the motion flows.

Fig. 1 .

 1 Fig. 1. The problem setup. The camera performs a planar motion on the ground plane. The motion parameters are the center of rotation C and the angle of rotation α. The images show a particular motion flow for a given point p0 in the image plane.

Fig. 2 .

 2 Fig. 2. Center estimation. c divides any line crossing it in two sets of points which are displaced by the motion flows in opposite halfspaces. Red lines shows the displacement vectors for points on the dashed segments. Left image shows a correct set of flows, whereas right image shows flows which suffer from drifting

Fig. 3 .

 3 Fig. 3. Angle estimation. Given an initial point p, by counting the number of steps required to cover the closed path, we are able to calculate the roto-translation angle

Fig. 4 .

 4 Fig. 4. Given a point p we detect the other intersection q of the two motion flows passing through p. This allows one to recover the two angles α1 and α2 which are used to rectify the ground plane. Points in blue represent the intersections of the two motion flows with line L.

Fig. 5 .

 5 Fig. 5. Estimated centers varying the b-spline parameters for two different planar motion. Red points represents inlier points, which are correctly distributed around the ground truth, while white points represents outliers. Notice that, in these experiments, the outliers are easily identified. Right image shows a close up of the two areas

Fig. 6 .Fig. 7 .

 67 Fig. 6. Angle estimation. The red curve represent the application of function Γᾱ(•, 1), starting from an initial point p. The segment connecting p to the estimated center is shown in white. It is used to count how many times the curve traverses it. Right image shows the estimation of the angle (in degree) w.r.t. the number of points used (i.e. the number of application of Γᾱ(•, 1))

Fig. 8 .

 8 Fig. 8. The function H1 and H2 for a particular experiment. Both functions are constituted by concentric flows around the centers of rotation. The interpolation allows one to detect values also for occluded region, such as the center of the image. Notice however that, at least for the function shown in blue, the interpolation fails near the external region of the image

Fig. 9 .

 9 Fig. 9. Left image shows the nodes used to recover the rectification constraints. Middle image shows the position of those nodes in the rectified image. Right image shows the rectified ground plane


Haruka Azuma
email: azuma@am.sanken.osaka-u.ac.jp

Yasuhiro Mukaigawa
email: mukaigaw@am.sanken.osaka-u.ac.jp

Yasushi Yagi
email: yagi@am.sanken.osaka-u.ac.jp

Spatio-Temporal Lifelog Using a Wearable Compound Omnidirectional Sensor

In this paper, we propose a spatio-temporal lifelog, that enables us intuitively to understand both the spatial and temporal aspects of a situation. We have designed a new wearable compound omnidirectional sensor comprising a hyperboloidal mirror and multiple paraboloidal mirrors. Using this sensor, visual information is recorded in all directions around the user. Objects are classified as foreground or background according to their distance. The stored log can be viewed using a browser that supports two different views, a spatial view and temporal view, which allow an effective visualization of the situation based on the classification. The user can switch between views depending on the purpose.

Introduction

Throughout human history, man has attempted to record his activities to enable these to be recalled at a later stage. Traditional recording methods, such as writing a diary or taking photos, are now maintained by digital devices. Recently there has been a trend to record all our daily activities as digital data without any human intervention. By storing such information automatically it would be possible to search for a person you once met or relive an ordinary day from the past as if you watch a movie of yourself.

Previous studies have focused on achieving an effective automatic lifelog. Microsoft Research developed a wearable camera containing other sensors such as a motion sensor [START_REF] Gemmell | Passive capture and ensuing issues for a personal lifetime store[END_REF]. The camera with a fish-eye lens covers almost all of the user's view. In one of the studies using the sensor, captured images are browsed with different sized icons ranked by frequency of appearance [START_REF] Smeaton | Content vs. context for multimedia semantics: The case of sensecam image structuring[END_REF]. Aizawa et al. [START_REF] Aizawa | Capture and efficient retrieval of life log[END_REF] used brainwave data to distinguish an image including the object which caught the user's attention. In all these studies, captured images have been used to re-create or complement a user's visual memory.

These previous lifelog systems have not, however, focused on spatial and temporal understanding of the event. Spatial information such as the directions of the surrounding people acts as a clue when recalling the past. Temporal information such as continuous relationships with different people is also useful if it is given all together. To develop a spatio-temporal lifelog, we need to acquire visual information that surpasses that which is detectable by the user's eye. Moreover, user foreground background Fig. 1. Model of the surrounding environment. The foreground contains a group of objects that could be elements included in an activity related to the user, while the background contains a group of objects specifying the location in which the user exists.

it is necessary to visualize such rich information in a way that allows the user to understand the situation intuitively.

In this paper, we propose a spatio-temporal lifelog using a compound omnidirectional sensor that can acquire visual information around a user and can classify the objects in the environment according to distance. We have designed a new wearable sensor that consists of a hyperboloidal mirror and multiple paraboloidal mirrors. We have also developed a browser with two different views aimed at enhancing spatial and temporal understanding, respectively.

Spatio-Temporal Lifelog

A spatio-temporal lifelog is a tool to aid the understanding of a situation from both spatial and temporal aspects. We need spatial information to perceive the situation surrounded by other objects such as friends or buildings. The directions from where we are to these objects and positional relations among the objects are key to recalling the situation precisely. Additionally to ascertain how situations have changed over time without spending too much time checking, we need to see all the changes in one glance. Temporal information specifies locations which change continuously as we walk and as a result of the motion of objects around us.

To clarify the requirements for a spatio-temporal lifelog, we first model the surrounding environment as illustrated in Fig. 1. The lifelog user is in the center, surrounded by two categories of objects, namely foreground and background objects, according to the context of association with the user's event. The fore-ground comprises a group of objects that could be elements included in an activity related to the user. On the other hand, the background comprises a group of objects specifying the location in which the user exists. In other words, objects are distinguished according to their roles with respect to the user. For example, while working in an office, a computer screen in front of you or a colleague sitting next to you would be considered part of the foreground. In contrast, the walls of the room surrounding you or some colleagues working farther from you would be background objects.

To build a spatio-temporal lifelog based on this environmental model, visual information from all directions around the user is necessary. The objects in the environment also need to be classified as foreground and background. Finally the acquired information needs to be presented in a meaningful way to enable spatial and temporal understanding of the situation.

Compound Omnidirectional Sensor

Sensor Design

When recording, it is necessary to capture the scene around the user together with information to classify objects as foreground or background. We distinguish foreground and background according to the distance from the sensor, with near objects being classified as foreground and far objects as background. This solution is valid where an important object such as a friend, who is talking to you and is immediately in front of you, is classified as foreground.

To obtain omnidirectional view and distance information, various catadioptric stereo systems have been proposed. Jang et al. [START_REF] Jang | Single-camera panoramic stereo system with singleviewpoint optics[END_REF] proposed a small sensor with two paraboloidal mirrors upside down, but which only has a single baseline for stereo matching. Several compound omnidirectional sensors have been proposed to allow multiple baselines for stable distance estimation. Kojima et al. [START_REF] Kojima | Calibration and performance evaluation of omnidirectional sensor with compound spherical mirrors[END_REF] developed a compound omnidirectional sensor with a spherical mirror in the middle and smaller ones around it. The problem with this sensor is that a complete perspective image cannot be produced because a spherical mirror does not maintain a single viewpoint. Ngo et al. [START_REF] Thanh | Robust and real-time rotation estimation of compound omnidirectional sensor[END_REF] proposed a sensor with seven paraboloidal mirrors, all of which have their respective single viewpoint. However, the telecentric lens used in this sensor resulted in the sensor being too big to be attached to a human body.

For use in a lifelog, we have designed a new wearable compound omnidirectional sensor, that consists of a hyperboloidal mirror, a concave paraboloidal mirror and four paraboloidal mirrors. The optical part of our sensor is a combination of two layers as shown in Fig. 2. The larger mirror, illustrated in (a), is hyperboloidal and is used to obtain a high resolution image. A ray going to a focal point of the hyperboloidal mirror, depicted as F h in (a), and which is a viewpoint, is reflected to a single point F h , which is another focal point. On the other hand, the smaller paraboloidal mirrors are used to supply multiple viewpoints and are illustrated in (b). A ray going to a focal point of the convex hyperboloidal mirror, depicted as F p in (b), is reflected vertically to the upper concave paraboloidal mirror, and then folded to the point F p , which is a focal point of the concave paraboloidal mirror. A folded catadioptric system [START_REF] Nayar | Folded catadioptric cameras[END_REF] adopted in layer (b) shrinks the overall sensor size.

The integrated optical schematic of the sensor is shown in Fig. 3. (a) is a crosssection view of the overall optical part and (b) is a view from the top of four paraboloidal mirrors, with diagonal lines indicating the cross section. A light ray L 1 from O(r, z) going to the focal point F h , which is also a viewpoint, is reflected to another focal point F h , and finally projected to P 1 on the image plane. On the other hand, a light ray L 2 going to the focal point F p is reflected in parallel to the optical axis z, then folded by the concave paraboloidal mirror, passes the focal point F h , which is also a focal point of the concave paraboloidal mirror, and finally reaches P 2 on the image plane. Multi-baseline stereo is achieved among the hyperboloidal mirrors and four convex paraboloidal mirrors.

Sensor Specification

Here we describe the setup of our compound omnidirectional sensor in detail. Figure 4 (a) is the overall view of the sensor. The height of the optical part is 50mm and the diameter of the central hyperboloidal mirror is 40mm. The overall height including the camera is 85mm. The weight of the sensor is only 50g because the mirrors are made of plastic with an evaporation coating. The camera has a 1/3 inch CCD with 480 by 640 pixels.

The vertical baseline between the hyperboloidal mirror and the convex paraboloidal mirrors is 30.2 mm long, while the shortest and longest horizontal baselines of the convex hyperboloidal mirrors are 8.3mm and 14.0mm, respectively. A vertical overlap of the field of view among the hyperboloidal mirror and convex paraboloidal mirrors is between -20.2 degrees and 11.5 degrees.

Figure 4 (b) shows an example of a captured image, with (c) depicting its panoramic image. The convex paraboloidal mirrors are visible at the four corners in the image in (b).

Near Object Detection

The objects can be classified by distance with the compound omnidirectional sensor. In principle, it is possible to calculate the precise depth by searching the corresponding points between mirrors [START_REF] Gluckman | Real-time omnidirectional and panoramic stereo[END_REF]. However we only need two classifications based on whether an object is near or far from the sensor. For this reason, we have adopted a method of discriminating between only two alternatives, near or far [START_REF] Kojima | Calibration and performance evaluation of omnidirectional sensor with compound spherical mirrors[END_REF].

If we assume that every ray to an image plane comes from infinity, there is no disparity among mirrors. The resolution of the camera defines a minimum distance that is sufficiently far and treated as infinity. Corresponding points between mirrors are calculated based on their positional relation, assuming that a ray from infinity reaches each mirror in parallel. Finally we know whether a ray comes from closer or further than the minimum distance by comparing the pixel values between these corresponding points. If the difference of the pixel values is less than a threshold, a far object is seen in the area. In contrast, if the difference is more than the threshold, there is a near object.

Figure 5 shows the correspondence between the hyperboloidal mirror and paraboloidal mirrors when every ray comes from infinity. The corresponding pixels are marked with the same color. Image size is normalized to fit the area of the mirrors. The inner zone of the hyperboloidal mirror corresponds to the outer Three examples of near object detection are shown in Fig. 6. The left column shows the inputs, while the right column illustrates the resulting perspective images where the area detected as near is marked in red. In the top image, where the person is at a distance of 4m, no near area is detected. In the middle and bottom images, where the person is at a distance of 2m and 50cm respectively, the area around the person is detected as near.

Lifelog Browser

Multiple Views

We now possess omnidirectional visual information around the user and a classification of foreground and background. Next, we focus on a means of browsing which is suitable for understanding the situation both spatially and temporally. The guiding principles we followed to build the browser are given below: • The changing situation in the foreground is recognizable at a glance without checking sequential images from a video. • The location where a person exists can be perceived.

-
Omnidirectional images are often viewed as a panoramic image, which achieves the same perspective view as that from our eyes. However, when we think of spatial understanding, a panoramic image does not suffice since it does not depict any spatial relations among objects in the image. To perceive a situation spatially, three-dimensional expression is needed. Meanwhile, to enable us to deal with a situation temporally without having to check a sequence of images, threedimensional expression with time as the third axis is crucial.

Consequently, we propose a browser with two different views that facilitate spatial and temporal understanding, respectively. The two views have been designed to enable the user to understand the situation intuitively from different aspects rather than presenting the whole omnidirectional image captured by the sensor. The user can switch between these views depending on the purpose. In addition, the viewpoint can be freely controlled to display the area the user wishes to see.

Spatial View

As shown in Fig. 7, spatial view depicts the background as a cylinder and foreground objects as billboards inside. (a) is our browser showing spatial view and (b) is an explanation thereof. The cylinder represents one of the backgrounds. Foreground objects appear and disappear with the automatic progression of time. The arrow at the bottom indicates the reference direction related to the position of the sensor. In this view temporal progress is merely shown as a series of images as in a movie. Dragging the cylinder with a mouse causes it to be rotated in any direction for better visibility. The buttons on the right of the window can be used to translate the viewpoint into a three dimensional space. In this view, the positional relations between the user and foreground or background objects, or the relations among foreground objects are easily recognized. By representing images as a three-dimensional form, the environment is visualized realistically and can be perceived intuitively. This view allows a user to concentrate on understanding the situation spatially, while some areas are designed not to be seen because of the shape of the displayed image.

Temporal View

The temporal view represents the background as a road and the foreground as billboards on it as shown in Fig. 8. (a) is our browser showing the temporal view and (b) is the explanation thereof. Temporal progress is expressed as depth, with later time on the far side. The background appears continuously in the road to inspire the user to recall the place where he was standing. Foregrounds appear as billboards walking along the background road. Multiple billboards represent a single person or object as they remained continuously around the user. In this view, as the user walks or as the foreground/background moves around the user whilst recording, a continuous change in textures is shown on the road or strings of billboards. As in the spatial view, the entire foreground and background can be rotated by dragging with a mouse, and the viewpoint can also be translated.

The main operation in this view is translation along the depth combined with temporal progress.

In this view, a user can easily recognize the changing situation at a glance without checking it in a movie. This view omits spatial representation in a threedimensional space to allow a user to focus on the temporal changes in the foreground and background.

Experimental Results

Using the compound omnidirectional sensor we recorded two situations in an indoor environment. For these experiments, the sensor was fixed to a tripod. Our compound omnidirectional sensor can classify objects as foreground and background, but the accuracy of the classification depends on the alignment. To demonstrate the effectiveness of the browser, foreground and background areas in the images are chosen manually. Foreground areas are demarcated with rectangles to be easily recognized. Figure 9 shows experimental results that demonstrate the advantage of the spatial view over panoramic images. The situation with a sensor in the center of a room and a person walking towards the sensor and then stepping away along the path indicated by the line, is illustrated in (a),. In the three panoramic images shown in (b), it is difficult to follow the person knowing his direction from the sensor. A further problem here is that his trajectory is segmented on the way. In contrast, the three representative images in the spatial view shown in (c) allow us to see the person passing on the left of the sensor. His movement is correctly recognizable in terms of the positional relation to both the sensor and the background.

Figure 10 shows the experimental results that demonstrate the advantage of the temporal view over panoramic images. The situation in which a sensor is located in a corridor and three people walk towards the sensor and then step away along the path indicated by the three different lines is illustrated in (a),. In the panoramic images shown in (b), we can see that three people appear around the sensor, but the three images are not enough to understand their motion. The representative images snapshot from the three different viewpoints in the temporal view are shown in (c). As observed in the first image, the continuous movement of multiple people is recognized at the same time with a single still image. In the second image we can see that one of the people is about to descend the stairs. The third image depicts a series of continuous motion of another person. As demonstrated in these examples, temporal view is useful to understand temporal change without time consuming viewing.

Conclusions

In this paper we have proposed a new spatio-temporal lifelog which enables us to understand a situation both spatially and temporally. We modeled the environment as a combination of foreground objects, which are directly related to a user's event, and backgrounds, which specify the location in which the user exists. We designed a new wearable compound omnidirectional sensor that can acquire visual information around the user and can classify the objects in the environment as foreground and background objects based on their distance. We developed a browser with both spatial and temporal views, which allows a user to understand a situation from two aspects.

Future work includes accurate calibration and alignment of the mirrors and the camera for better accuracy of near and far object classification. Focusing will be improved by replacing the current lens with another which has the wider focal depth. We also plan to add other views to the browser, making use of omnidirectional visual information and separation of foregrounds and backgrounds.

This research was partially supported by the Ministry of Education, Science, Sports and Culture, Grant-in-Aid for Scientific Research (S).

FhFig. 2 .Fig. 3 .Fig. 4 .

 234 Fig. 2. Configuration of two layers: (a) the layer with a hyperboloidal mirror, and (b) the layer with paraboloidal mirrors. Both configurations retain a single viewpoint (F h /Fp) and converge a ray to a single point (F h /F p ).

Fig. 5 .

 5 Fig. 5. Correspondence between the hyperboloidal mirror and surrounding paraboloidal mirrors when every ray is from infinity.

Fig. 6 .

 6 Fig. 6. Three examples of near object detection. The left column shows the inputs, while the right column illustrates the resulting perspective images where the area detected as near is marked in red.

Fig. 7 .

 7 Fig. 7. An example of spatial view: (a) the browser, and (b) a description thereof. A user can easily recognize the local relation between a person and the foreground/background objects or among foreground objects.

Fig. 8 .

 8 Fig. 8. An example of temporal view: (a) our browser, and (b) the description thereof. A user can easily recognize the changing situation at a glance.

Fig. 9 .

 9 Fig. 9. Experimental results of spatial view: (a) the situation, (b) the panoramic images and (c) the spatial view. The compound omnidirectional sensor is in the center, a person approaches the sensor and then steps away along the line.

Fig. 10 .

 10 Fig. 10. Experimental results of temporal view: (a) the situation, (b) panoramic images, and (c) the temporal view. The compound omnidirectional sensor is in the center, and three people approach the sensor and step away along the lines.