Destination Prediction by Trajectory Distribution Based Model

Philippe C. Besse, Brendan Guillouet, Jean-Michel Loubes, and Franc ¸ois Royer

Abstract-In this paper we propose a new method to predict the final destination of vehicle trips based on their initial partial trajectories. We first review how we obtained clustering of trajectories that describes user behaviour. Then, we explain how we model main traffic flow patterns by a mixture of 2d Gaussian distributions. This yielded a density based clustering of locations, which produces a data driven grid of similar points within each pattern. We present how this model can be used to predict the final destination of a new trajectory based on their first locations using a two step procedure: We first assign the new trajectory to the clusters it mot likely belongs. Secondly, we use characteristics from trajectories inside these clusters to predict the final destination. Finally, we present experimental results of our methods for classification of trajectories and final destination prediction on datasets of timestamped GPS-Location of taxi trips. We test our methods on two different datasets, to assess the capacity of our method to adapt automatically to different subsets.

Index Terms-Trajectory Classification, Final Destination Prediction.

I. INTRODUCTION

M ONITORING and predicting road traffic is of great importance for traffic managers. With the increase of mobile sensors, such as GPS devices and smartphones, much information is at hand to understand urban traffic. In the last few years, a large amount of research has been conducted in order to use this data to model and analyze road traffic conditions. The aim of this paper is to tackle the issue of predicting the destination of vehicles given a prefix of their trajectory. This problem has been the subject of a Kaggle challenge entitled "ECML/PKDD 15: Taxi Trajectory Prediction (I)" [START_REF]KAGGLE data set ecml/pkdd 15: Taxi trajectory prediction (1)[END_REF].

The observations are time-stamped locations that correspond to the different positions of vehicles moving within a city monitored at different observation times. When dealing with a dataset composed of trajectories, the difficulty lies in the fact that the data convey both spatial information (locations of the vehicles on the map of the city) and temporal information (for each vehicle, the locations are indexed by time, which creates a sequence of locations that compose a full trajectory). Hence the data have a spatio-temporal structure that must be taken into account in order to model their evolution while the trajectories of the destination points to be predicted are unknown. Vehicle trajectories are also constrained to a road network which makes their time progression very irregular. Locations of vehicles can be seen as two-dimensional data in R 2 , that have to be compared to one another, taking into account characteristics of the trajectories they belong to, such as origin and destination.

In this paper, we propose a method that relies on a distribution based model for the trajectories. We first focus on the temporal structure of the data and gather the locations into clusters of points that belong to similar trajectories. For this, we rely on a distance that takes into account the geometric properties of trajectories which was developed in a previous work [START_REF] Besse | Review and perspective for distance based trajectory clustering[END_REF]. Then, we model the observations within each obtained cluster by the realisation of a random variable that must be estimated with a distribution on R 2 . This estimation step is achieved by considering a mixture of 2-dimensional Gaussian distributions fitted to the data by a maximum likelihood procedure. Thus, each cluster of trajectories corresponds to a parametric distribution model obtained by a mixture of Gaussian distributions. Using this procedure, we obtain a distribution model for points based on the assumption that they belong to a cluster of trajectories. Forecasting the destination of vehicles is a two step procedure :

• we first attribute the observed path of these vehicles to a cluster of trajectories and • then extract from the trajectories within the cluster a feature that stands for the final destination point. Hence using the learning set, we obtain a density classification method based on preliminary trajectories clustering. In addition to the forecast properties of this model that will be analysed on the taxis data set, this methodology provides a probabilistic model for spatio-temporal analysis of vehicle flows. It enables the extraction of distribution mobility patterns in an vehicle transportation system that in order to understand urban mobility and flow.

The paper is organised into the following sections. Section II is devoted to the presentation of the data and the related work. In Section III we present how we obtained trajectories clustering by using a proper distance that takes into account the spatial properties of vehicle trajectories. Section IV describes the issue of clustering points into clusters of trajectories with a mixture of Gaussian distributions and how to use these models to predict the final destination of taxi trip. Finally, in Section V, we present the experimental results and the performance of our models.

II. PRESENTATION OF THE FORECAST PROBLEM AND RELATED WORK

Consider vehicles' location data that consist of locations p(t) ∈ R 2 observed at several observation times t that may differ for each vehicle's path. As an example of such data, through the paper we will test our procedure using two different datasets. The first contains over 11 million taxi-GPS samples of approximately 500 taxis collected over 30 days in the San Francisco, United-States [START_REF] Piorkowski | CRAW-DAD data set epfl/mobility (v. 2009-02-24)[END_REF]. The second contains more than 83 millions taxis-GPS data points describing July 2013 to June 2014 for all 442 taxis in circulation in Porto, Portugal. This dataset has been provided for a Kaggle Competition [START_REF]KAGGLE data set ecml/pkdd 15: Taxi trajectory prediction (1)[END_REF]. This dataset is also composed of metadata associated to the taxi trips, such that client, taxi stand or taxi identification but the dataset in San Francisco only got the taxi identification. Hence, we deliberately do not use these attributes because we want our method to be easily adapted from one dataset to another. We only use locations in R 2 and the associated timestamp.

Our goal is to be able to determine the final destination point by observing the beginning of a trajectory. This issue is common to various area and data, animal migrations [START_REF] Gaffney | Trajectory clustering with mixtures of regression models[END_REF], VideoFrames applied on Robotic [START_REF] Kruse | Acquisition of statistical motion patterns in dynamic environments and their application to mobile robot motion planning[END_REF], Human [START_REF] Vasquez | Motion prediction for moving objects: a statistical approach[END_REF] Vehicle on crossroads [START_REF] Hu | A system for learning statistical motion patterns[END_REF], or GPS data to study behaviour in Urban Commercial Complexes [START_REF] Kim | Adnext: a visit-pattern-aware mobile advertising system for urban commercial complexes[END_REF]. But vehicle trajectories are very different objects, they are constrained to a road network and have very irregular time progressions. The study of destination prediction requires comparing the information of previous trajectories with the current location of trajectories in order to identify the destination. Many authors have already discussed this issue. In the winning solution of the Kaggle-ECML/PKDD discovery challenge on taxi destination prediction, De Brebrisson et al. [START_REF] De Brébisson | Artificial neural networks applied to taxi destination prediction[END_REF] used a multi-layer perceptrons neural network on features vector composed of coordinates of beginnings of trajectories, and diverse context information, such as the departure time, the driver id and client information. Their training set has been built to match the trajectories in the test set's competition. If this solution can easily be adapted to other datasets, it implies a new training for some tested trajectories using more location information. Moreover, neural network scores are hard to interpret and can not be used to better understand the characteristics of the dataset. Krumm et al. [START_REF] Krumm | Predestination: Inferring destinations from partial trajectories[END_REF] and Ziebart et al. [START_REF] Ziebart | Navigate like a cabbie: Probabilistic reasoning from observed context-aware behavior[END_REF] also used external information, in addition to historical trajectories, such as travel time, trajectory length, accident reports, road condition, and driving habits. They incorporate this information into Bayesian inference to compute the probabilities of predicted destinations. Parteson GRAPH et al. [START_REF] Patterson | Inferring highlevel behavior from low-level sensors[END_REF] also used a Bayesian method to predict destination but for specific individuals based on their historical transport modes. The main idea of these studies is to use the external information to enhance the quality of the prediction. It then becomes dependent on the presence of this information and is inapplicable without them.

Monreale et al [START_REF] Monreale | Wherenext: a location predictor on trajectory pattern mining[END_REF] built a decision tree, named T-pattern Tree based on extracted movement patterns and predicted the next location of a new trajectory finding the best matching path in the tree. Tiesyte and Jensen [START_REF] Tiesyte | Similarity-based prediction of travel times for vehicles traveling on known routes[END_REF] proposed a nearestneighbour trajectory method that utilised distance measures to identify the historical trajectory most similar to the current partial trajectory.

Finally, most of the work dealing with destination prediction issue uses probabilistic methods based on the location to identify the most probable location after creating the probability model. Among them, the Markov model has been widely applied in predicting destinations. Ashbrook et al. [START_REF] Ashbrook | Using gps to learn significant locations and predict movement across multiple users[END_REF] find potential destinations by clustering GPS data, then predicting destinations from these candidates based on Markov models trained to find the next most likely destination based on those recently visited. Gambs et al. [START_REF] Gambs | Next place prediction using mobility markov chains[END_REF] determined the destination by using the mobility Markov chains from the sequence of the POI (Point of Interest) to create the model. Simmon et al. [START_REF] Simmons | Learning to predict driver route and destination intent[END_REF] also built the probabilistic model through observation of the drivers' habits. All these studies build prediction based on habits of one or a group of specific individuals based on their historical trips. But they require knowing the identity of the driver. Xue et al. [START_REF] Xue | Destination prediction by sub-trajectory synthesis and privacy protection against such prediction[END_REF] proposed a method which decomposes all available patterns into subtrajectories of neighboring locations. The subtrajectories are assembled into synthesized trajectories. Then, they build the Markov model, which quantifies the correlation between adjacent locations. The main drawback of both bayesian inference and the Markov model is in establishing how well they discretize the space. Either they use the true road network, [START_REF] Patterson | Inferring highlevel behavior from low-level sensors[END_REF], [START_REF] Ashbrook | Using gps to learn significant locations and predict movement across multiple users[END_REF], [START_REF] Simmons | Learning to predict driver route and destination intent[END_REF], which requires significant amount of extra work to map the GPS data to the graph network, or they use a grid of square cells [START_REF] Krumm | Predestination: Inferring destinations from partial trajectories[END_REF], [START_REF] Ziebart | Navigate like a cabbie: Probabilistic reasoning from observed context-aware behavior[END_REF], [START_REF] Xue | Destination prediction by sub-trajectory synthesis and privacy protection against such prediction[END_REF] which is a rough representation of the space and produces results dependent on the choice of the discretization of the grid.

Choi and Hebert [START_REF] Choi | Learning and predicting moving object trajectory: a piecewise trajectory segment approach[END_REF] present a Markov model based on segments of the trajectories, where latent segments are obtained by clustering segments of past trajectories. New trajectories are then modelled as a concatenation of segments, which are assumed as noisy realisations of the latent segment. One of the major drawbacks of this method is that it is used for short term prediction (at most 10 seconds ahead). Wiest el al. [START_REF] Wiest | Probabilistic trajectory prediction with gaussian mixture models[END_REF] proposed a probabilistic trajectory prediction based on two types of mixture models. They also predict the vehicles trajectories only several seconds into the future.

III. MODEL TO CLUSTER TRAJECTORIES

In this section we describe how we cluster trajectories. We first recall the definition of trajectory used in this paper.

Definition 1. A trajectory T i is defined as T i : [(p i 1 , t i 1 ), . . . , (p i n i , t i n i )], where p i k ∈ R 2 , t k ∈ R ∀k ∈ [1, . . . , n i ], ∀n i ∈ N and n i is the length of the trajectory T i .
To cope with the sampling issues of trajectories, we first complete, when required, the locations between p i j (at time t i j ) and p i j+1 (at time t i j+1 ) by the piece wise linear representation between each successive location p i j and p i j+1 resulting in a line segment s i j between these two points. This new representation is called the piece wise linear trajectory. In this representation, no assumption is made about time indexing of segment s i j . Definition 2. A piece wise linear trajectory is defined as

T i pl : (s i 1 , . . . , s i n i -1 )
, where

s i j = [p i j , p i j+1 ] ∈ R 4
and n i pl is the length of the piece wise linear trajectory.

The length of the PL-trajectory n i pl is the sum of the lengths of all segments that compose it :

n i pl = j∈[1...n i -1] p i j p i j+1 2 .

TABLE I NOTATION

T The set of trajectories

T i
The i th trajectory of set T

n i Number of locations in trajectory T i p i j
The j th location of

T i t i j
The time index of location p i j l i Label of the i th Trajectory

T m = {T i |l i = m} Set of trajectories in cluster m P m = {p i j |l i = m} Set of points in cluster m C(T ) = {T 1 , . . . , T K } Set of cluster of trajectories K Number of clusters of trajectories in C(T ) l i j Label of the location p i j P m n = {p i j |l i j = m, l i = n} Set of points in cluster (m, n) C(P m ) = {P m 1 , . . . , P m k m } Set of points cluster k m Number of clusters of points in C(P m ) #S Number of elements in the set S
In a previous paper [START_REF] Besse | Review and perspective for distance based trajectory clustering[END_REF], we proposed a method to cluster trajectories based on the behaviours of the users. This clustering is obtained by hierarchical clustering with the ward linkage criterion based on a distance between trajectories, the Symetrized Segment-path-Distance.

Definition 3. Symmetrized Segment-Path Distance D SSP D (T 1 , T 2 ) = D SP D (T 1 , T 2 ) + D SP D (T 2 , T 1 ) 2 .
where

D SP D (T 1 , T 2 ) = 1 n 1 n1 i1=1 D pt (p 1 i1 , T 2 ).
and where D pt is the distance from a point to a trajectory as defined below.

Definition 4. P oint -to -Segment distance.

D ps (p 1 i1 , s 2 i2 ) = p 1 i 1 p 1proj i 1 2 if p 1proj i 1 ∈ s 2 i 2 , min( p 1 i 1 p 2 i 2 2, p 1 i 1 p 2 i 2 +1 2)
otherwise.

Where p 1proj i1 is the orthogonal projection of p 1 i1 on the segment s 2 i2 .

This distance compares trajectories as a whole, regardless of their time indexing or the number of locations that compose them. It enables us to produce clusters of trajectories describing the traffic flow of the trajectory set T . The partition, C(T ) of T shows the K main paths taken by the users. The number of clusters K, is not fixed. It depends on the dataset, or the precision we want to use in order to described it. In Section V, we discuss the choice of K for both datasets according to the values of the different quality criteria described in the same Section.

IV. A PROBABILISTIC MODEL FOR TRAJECTORY CLASSIFICATION

After obtaining clusters of trajectories that discriminate the main patterns of the traffic flow in the city, we aim to predict the final destination of a vehicle for which we only observe the beginning of its path. Hence, We observe a succession of locations in R 2 . To assign these points to a cluster of trajectories, we model the clusters by a mixture of 2d Gaussian distributions. Thus, for each cluster, we obtain a Gaussian likelihood estimated using only the data belonging to this cluster. The observed locations will then be assigned to the most likely cluster according to these different likelihoods.

A. Points Partitioning within Clusters of Trajectories

Recall that we obtain a set of clusters of trajectories, C(T ).

For a new trajectory, we want to be able to assign it to the cluster it most likely belongs to. For this purpose, we build a Gaussian mixture model for every cluster of trajectories T m ∈ C(T ) from the set of all the points P m that compose these trajectories. A Gaussian mixture model assumes that all points from P m are generated from the sum of k m Gaussian distributions φ , which are, in our case, 2-variate Gaussian distributions.

Definition 5. A Gaussian Mixture Model is a weighted sum of k m component Gaussian densities as given by the equation,

Φ m (p) = Φ(p|Θ m ) = k m k=1 ω m k • φ m k (p),
where ω m k is the mixture weight, i.e. the prior probability for any point p belonging to the k th cluster, such that Each component density is a Gaussian density on R 2 . Definition 6. The density function of a normal distribution, φ m k , is defined as,

φ m k (p) = φ(p|µ m k , Σ m k ) = 1 (2π) 2 |Σ m k | e -1 2 (p-µ m k ) tr •(Σ m k ) -1 •(p-µ m k ) ,
where µ m k ∈ R 2 and Σ m k ∈ R 2x2 are respectively the location and the covariance matrix of φ m k , and

|Σ m k | is the determinant of the matrix Σ m k . Θ m = {ω m 1 , µ m 1 , Σ m 1 , . . . , ω m k m , µ m k m , Σ m k m } is the list of parameters of the GMM distribution Φ m .
To evaluate the parameters Θ m , we use the maximum likelihood estimation. Its aim is to find the parameters which maximise the likelihood function of Φ m , given the training set P m . The GMM likelihood, can be defined as,

L(Θ m |P m ) = p∈P m Φ m (p). (1)
The maximum likelihood estimators, Θ m M L , are the parameters which maximise the GMM likelihood function.

Θ m M L = arg max Θ L(Θ m |P m ). (2)
For each GMM, the number of components k m is set to the value which maximise the criterion information BIC

= -2 ln L(Θ m |P m ) + k ln(|P m |). k m = arg max k BIC(k). (3)
The complete set of trajectories is then modelled by the set of K GMM's, one for each set of points, P m . Each of these sets has been partitioned into k m groups : C(P m ) = {P m 1 , . . . , P m k m }. Using this modelling procedures, we obtain several cluster of locations, each one corresponding to a mode of the estimated Gaussian mixture distribution. We got a density based clustering of the cloud points, which produces a data driven grid of similar points within a cluster of trajectories. Now that we have described the space, we want to use the model to predict the final destination of new trajectories. For that, we want to be able to assign the new trajectory to the cluster it most resembles.

B. Classification of Trajectories

Here, We present our method for classification of trajectories. For a new trajectory T c , we want to assign it to the cluster of trajectories it most likely belongs. For this purpose we compute the simple score, s m (T c ) for all the GMMs Φ m . The score is the value of the likelihood function of Φ m given the points that compose the trajectory T c . It represents how likely the trajectory T c belongs to the cluster m. Definition 7. The simple score, s m (T c ), for a trajectory, T c , to be assigned to the cluster m is defined as:

s m (T c ) = L(Θ m M L |T c ) = P (T c |Θ m M L ) = p c j ∈T c Φ m (p c j |Θ m M L )
In this way, we can assigned the trajectory to the cluster with the highest affinity score.

l c guess = max m∈[1...K] s m (T c ). (4)
We highlight the fact that this method enables us to compute a score for the trajectory and for each cluster. The trajectory is not attributed to one cluster and one cluster only. This is relevant because when only a few points of the trajectory are known, we cannot always be totally certain of the final destination. Several destinations are possible. Hence the score computed is a probability that the trajectory belongs to the cluster of trajectories.

C. Complete Model

We want to test the influence of auxiliary variables on the quality of our classification method, such hour of the day, or day of the week, during which the trip takes place. The likelihood score as defined in Definition 7 does not take into account contextual information. However, we can assume that prior knowledge may help to discriminate the trajectories. Indeed, a path may more likely be taken than an other at a given hour of the day or day of the week. We look forward to verifying this hypothesis by including auxiliary weights. For this we define a new complete score taking into account the following weights. Definition 8. The complete score, s m c (T c ), for a trajectory, T c , to be assigned to the cluster m is defined as:

s m c (T c ) = L ap (Θ m M L |T c , E c ) = P (E c |Θ m M L )P (T c |Θ m M L ) = α(m, h c , d c ) p c j ∈T c Φ m (p c j |Θ m M L )
The definition of the complete score( 8) is generic. h c ∈ [0, . . . , 23] and w c ∈ [1, . . . , 7] are respectively the hour of the day and the day of the week at which the trajectory T c begins. This information can be interpreted in different manners which will result in a different value for the auxiliary weight α(m, h c , w c ) according to the information we are taking into account. We define three different weights:

• The Empiric weight describes the distribution information of the trajectory cluster.

a emp (m) = #T m #T .
• The Weekday weight describes the distribution information of the trajectory cluster at a given day of the week,

a wd (d, m) = #{T c | T c ∈ T m , d c = d} #{T c | T c ∈ T , d c = d} .
• The Hours weight describes the distribution information of the trajectory cluster at a given hour of the day,

a h (h, m) = #{T c | T c ∈ T m , h c = h} #{T c | T c ∈ T , h c = h} .
The auxiliary weight α(m, h c , w c ) is the product of any combination of these weights.

D. Model for Final Destination Prediction

We present here how our model can be used to predict the final destination of the user trips. We have defined, section IV-B, a simple score and a complete score for each trajectory to belong to a cluster of trajectory. Hence we can assign the new trajetory to the clusters it most likely belongs. We can then use the information from the trajectories that compose these clusters to predict the final destination of the new trajectory. From this, we define two different methods for predicting final destination.

On the one hand, we consider only the trajectories from the cluster m with the highest score.

d pred1 (T c ) = 1 #T m i s.t. l i =m p i n i = d m , s.t. m = l c guess , (5)
where d m is the mean of the locations of all final destinations of the trajectories in cluster T m . On the other hand, in order to take advantage of the fact that the trajectory T c is not strictly assigned to one and only one cluster, we define, d pred2 , as a weighted sum of the mean final destination of every cluster.

d pred (T c ) = K m=1 s m (T c ) K k=1 s k (T c ) • d m = K m=1 s m w (T c ) • d m , (6)
where s m w (T c ) is the weighted affinity score of s m (T c ) The complete score can be used instead of the simple score to take into account the effect of the auxiliary variable in the final destination prediction.

V. EXPERIMENTAL RESULTS

In this section we present experimental results to evaluate both classification and final destination's prediction methods. To evaluate prediction error, we use the Haversine Distance (see Definition 11), which is the evaluation metric used in the Kaggle competition [START_REF]KAGGLE data set ecml/pkdd 15: Taxi trajectory prediction (1)[END_REF]. The Haversine Distance measures distances between two points on a sphere based on their latitude and longitude. We use a 10-cross validation method to calculate this error by learning on 90% the data: the training set T train , and forecasting the remaining 10%: the test set T test . The error forecast is the average for all the training sets. We repeat this operation ten times, such that every set has been considered as the test set, to ensure a more accurate estimation of model prediction performance.

To evaluate our method during trajectory completion, we introduce the definition of a partial trajectory, below. A ptrajectory, T i (p) of a trajectory T i is a subset of this trajectory such that the length of the piecewise representation of T i (p) is at most p times the size of the length of the piecewise representation of T i . Definition 9. The p-trajectory T i (p), ∀p ∈ [0, . . . 1] is defined as the trajectory:

T i (p) = ((p i 1 , t i 1 ), . . . , (p i n i (p) , t i n i (p) )) s.t. n i pl (p) n i pl ≤ p
where n i (p) is the number of locations that compose the ptrajectory T i (p).

A. Data And Clustering Results

To analyse our result and test its scalability, we test our model on two different subsets. The first one is a subset of taxi trajectories from San Francisco [START_REF] Piorkowski | CRAW-DAD data set epfl/mobility (v. 2009-02-24)[END_REF]. It is composed of 4.127 trajectories, all starting from the Caltrain Station and ending in an area of size 6, 327 × 6, 827 km in the center of the city. The second is composed of 19.423 trajectories from taxis in the center of Porto [START_REF]KAGGLE data set ecml/pkdd 15: Taxi trajectory prediction (1)[END_REF], leaving from the Sao Bento Station and ending in a delimited area of size 8, 116 × 8, 068 km. These two datasets are displayed on the left in Figure 1 and Figure 2.

These two datasets are quite different. In San-Francisco, the road network looks like a grid, most of the streets are either parallel or perpendicular to each other. In Porto the network is more irregular. These differences will enable us to test the scalability of our method and its capacity to be adapted on different datasets.

On figure 1 and2 we can see the results of the clustering described in section III on the data sets from San Francisco and Porto. We can see that the clustering obtained restuls in a group of trajectories tracing the same path from the selected departure points.

B. Trajectory Classification

To evaluate the quality of our classification, we observe the percentage of trajectories that have been assigned to their true cluster.

Definition 10. The quality criterion, Q class , for the classification, is the percentage of well classified p-trajectories , ∀p ∈ [0, . . . 1], defined as:

Q class (p) = #{T i (p)|l i guess = l i , T i ∈ T test } #{T i (p)|T i ∈ T test }
In Figure 3, we can observe the percentage of well classified trajectories, Q class , for p = 1, i.e, when the trips are completed, with respect to the number of clusters of trajectories. Its value does not necessary decrease when the number of clusters increases. When the number of clusters of trajectories is low, the points that compose the trajectory are scattered. Hence the clusters of points found with the Gaussian Mixture have covariance matrices with high value resulting in low likelihoods and low scores. This means that a trajectory can have a low score with respect to its correct cluster. This explains why increasing the number of clusters of trajectories does not always decrease the quality criterion Q class . Hence, it can be used as a good criterion to know how well the number of clusters chosen describes the traffic flow of the studied dataset. For the trajectories in San-Francisco, the highest value is 5 clusters, but the second highest value is obtained for 20 clusters. For the trajectories in Porto, the two highest values are found for 20 and 35 clusters.

We can observe that the percentage of well classified trajectories is always higher than 85% for trajectories in San Francisco and higher than 91% for trajectories in Porto even for 100 clusters. However, we do not achieve more than 95% of correct classification for both datasets. But our method does not strictly assign one cluster to each trajectory, but a score to each trajectory. We can see that the percentage of correctly classified trajectories significantly increased when we looked at the best-2 and best-3 predictions. For trajectories in Porto, this score is almost always greater than 99% for each number of clusters.

Hence, It is relevant to look at the score of the correct cluster, and not only the classification rate. For this purpose, we consider the ROC(Receiver operating characteristic) curves of a one-vs-all classification for every cluster that constitutes a clustering result.

In Figure 4, we can observe these ROC curves and their AUC (Area Under Curves) for the clustering fixed to 25 clusters for San Francisco, and 45 for Porto. All AUC are greater than 0.90, and 17 are greater than 0.95 for San Francisco. All but three are greater than 0.90 and 28 are greater than 0.95 for Porto. These results show that even if some trajectories are not assigned to its correct clusters, the simple scores for their correct cluster is always elevated.

C. Final Destination Prediction

In this section, we present the results of our method for the prediction of taxi trips destination. To evaluate our method we used the mean of the Haversine Distance, which measures distances between two points on Earth based on their latitude and longitude.

Definition 11. The Haversine Distance, D H between two locations d1, d2 ∈ R 2 is defined as:

D H (d1, d2) = 2 • r • arctan a 1 -a a = sin 2 y 2 -y 1 2 + cos(y 1 ) cos(y 2 ) sin 2 x 2 -x 1 2
where (x 1 , y 1 ) and (x 2 , y 2 ) are the longitude and latitude of d1 and d2 respectively and R = 6371(km) is the radius of the Earth. Hence, the Haversine distance returns the distance in Km between two locations on Earth.

We can then define the quality criterion for the prediction of the final destination as the mean of the Haversine distance between the true location of the final destination, p c nc , of the trajectory, T c , and the location of the prediction, d pred (T c ).

Definition 12. The quality criterion, Q pred , is defined as :

Q pred (p) = T c ∈T D H (d pred (T c ) -p c nc ) #T
In figure 5 we can observe the results of the quality criterion, Q pred , according to the trajectory completion. We compare the results of the two prediction methods, pred 1 and pred 2 . Both datasets are displayed, San Francisco(a), on the left and Porto(b), on the right, for a number of clusters of 25 and 45 respectively. For San-Francisco, we can observe that the second method gives best results especially at the beginning of the trajectories where Q pred is 400 meters better using pred 2 . As the trajectories progress, the results continue to be better with pred 2 but the difference between the two methods decreases and after 50% of trajectory completion, the difference is less than 50 meters. This is expected because the more locations we know for a trajectory, the more confidently we can assign the trajectory to one cluster and one cluster only. Hence, the more locations we know for a trajectory, the more closely the results are using the two methods. For trajectories in Porto, if Q pred also gives better results with pred 2 , the difference between the method is insignificant. Nevertheless, we will still use the pred 2 to compare results according to the number of clusters.

In Figure 6, we can look at the same quality criterion, Q pred , according to trajectory completion. We display these results for different numbers of clusters from 0 to 100. For dataset in San Francsico (a), and for trajectory completion between 0% and 50% the bests results are found for 5 and 10 clusters. At these completion rates, the trajectories are more easily assigned to their correct cluster, leading to best results. For completion rates between 50% and 100% the results found with 5 and 10 clusters are the worst, because these numbers of clusters do not well enough describe the space. The same conclusion can be made for 15 and 20 clusters. UP Until 70% of trajectory completion, there is no strong differences for a number of clusters between 25 and 100. When all the trajectories are completed, the more clusters we have, the more precise the prediction is. However, we have a gain of precision of only 200 meters between 25 and 100 clusters. The more clusters of trajectories we have, the more Gaussian Mixture we need to estimate. Hence 25 clusters of trajectories is the best compromise to well describe the space and to do so within a reasonable computation time. For Porto dataset(b), the worst results are found for a number of 5 clusters, for trajectory completions between 0% and 15% and between 40% and 100%. For trajectory completions from 35% to 55%, the best results are for a number of clusters of 10 and 15, but they yield bad results after 65% and 80% trajectory completions. The same conclusion can be made for a number of clusters between 20 and 40. We can observe that the results stabilise when the number of clusters increases from a number of cluster of 45. The difference of Q pred value for a number of clusters between 45 and 100 does not exceed 40 meters for trajectory completions from 0% to 80%. Similarly to San-Francsico, when all the trajectories are completed, the more clusters we have, the more precise the prediction is, but the gain of precision is low, only 100 meters between 45 and 100 clusters. Hence, 45 is the best choice for the Porto dataset.

D. Effect On Classification and Prediction

In Figure 7 and Figure 8, we observe the effect of different auxiliary weights described in Section IV-C on both the quality criteria Q class and Q pred . We display the differences of these criteria with the different weights and the same criteria with no weights according to trajectory completion. For San-Francisco, we display the results for 25 clusters. We can observe that all weights result in an improvement on both the quality of the classification and the prediction of the final destination in the first part of the trajectories, for trajectory completion between 0% and 35% -45%. The mix of all the weights is yields the best results. The improvement of classification continues until 8% when trajectory completion is at 25% and the improvement of the Q pred criteria is 225 meters when the trajectory starts and 100 meters at 25% of trajectory completion. The curves of both quality criteria, Q class and Q pred are not perfectly correlated. This is expected because Q class shows the rate of correct classification, while Q pred displays the prediction quality found with pred 2 , which uses information of different clusters and not only the first predicted cluster for the prediction. Beyond these completion rates, the auxiliary weights deteriorate according to the different quality criteria values. This means that when we have little information about the location of the trajectory, context information help to improve the destination prediction. Whereas when we have sufficient information about the trajectory location, we can confidently predict the correct clusters of trajectories the new trajectory most likely belongs to. Hence, adding auxiliary weight information deteriorates the result.

The results obtained are different using the Porto dataset. The results for Porto datasets are displayed for 45 clusters. The Fig. 9. Exemple of final destination prediction for a taxi trip In San Francisco different weights improve the prediction, and the mix of all weights yields the best results, but the improvement is always less than 30 meters which is much less significant than with the San-Francisco dataset. Similarly, the classification is never improved more than 3%. Taxi trips in Porto are less influenced by auxiliary variables than taxi trips in San Francisco.

In Conclusion, we have seen that our method gives similar results in trajectory classification, Section V-B, and in prediction of final destination, Section V-C, for studied datasets of trajectories in San Francisco and Porto. Taking into account the differences between the structure of the road network of theses two cities proves that our method can be adapted to different datasets, without requiring prior study of the dataset. However, the effect of auxiliary variables is different from one dataset to one another. These results show that the behaviour of the drivers differs from one city to another. It could help traffic managers to better understand the traffic flow of a city.

We have tested our methods with the test dataset from the Kaggle challenge. This competition has completed but we can still submit an entry to see our score. Our final results produce a mean error of 2, 82021 kilometres, while the best result was 2, 03489 kilometres. Hence we have a difference of 800 metres relative to the winning solution. This is a promising result because our model has been trained on a sample of trajectories which doest not strictly match the trajectories in the test set. Hence our model can be re-used directly for a different test dataset, and can also be used to predict the destination within the same trajectory, without requiring a new training. In Figure 9, we can observe an example of how our model for final destination prediction works. The probability of different possible destination points for a trajectory at 6 different percentages of its trip accomplishment is displayed in this figure. At each moment, the cluster of trajectories with its corresponding final destination and simple score are displayed. The more likely the trajectory belongs to a cluster, the more visibly this cluster is displayed on the plot.

VI. CONCLUSION

In this paper, we proposed a data-driven method to predict the final destination of vehicle trips using a statistic learning procedure. Vehicle trajectories differ from other trajectories in that they are constrained to a road network, which differs from one place to another, and directly influences the behaviour of the users. The learning step of our method follows a twostep procedure which enables to capture the behaviour of the user. It first models the main paths taken by the users by clustering their complete trajectories. Then, it models main traffic flow patterns within each trajectory's cluster by a mixture of 2d-Gaussian distributions. This yields a data driven grid of locations which describes the all space. This model is finally used to predict the final destination of vehicle trips, by assigning the trajectory to the path to whom it belongs the most likely and extracting information from trajectories who follow this path. This prediction is based on the initial location of the trajectory. Since we model the whole path, the prediction can be accomplished at any time during trajectory completion. Such method is applied on two different datasets: trajectories of taxi trip moving on two different road networks from San-Francisco, United-State and from Porto, Portugal and proves that such predictions based on the structures of the paths, compete with methods very complex and not easily tractable such as deep learning methods. Hence we propose a new description of road traffic that can be used for other research. For example, we can use different information from trajectories inside the clusters to short term prediction, estimate arrival time, or detect abnormal behaviour. Our model provide a better understanding of behaviours of the cars drivers by pointing out the main paths. Hence it can help organise trip distribution of a city.

k m k=1 ω m k = 1 ,

 1 and φ m k (p), i = 1, . . . , k m are the component Gaussian densities.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Caltrain Station, San Francisco Dataset and its Partitioning in 25 clusters

Fig. 3 .

 3 Fig. 3. Percentage of Trajectories Correctly Classified According to Number of Clusters. Compare Best-3 Prediction.

Fig. 4 .

 4 Fig. 4. Roc Curves and AUC for every clusters of trajectories.

Fig. 5 .

 5 Fig. 5. Mean Error of Final Destination Prediction According to Trajectory Completion. Compare Method 1 and 2.

Fig. 6 .

 6 Fig. 6. Mean Error of Final Destination Prediction According to Trajectory Completion. Compare Number of Cluster for Method 2

Fig. 7 .Fig. 8 .

 78 Fig. 7. Improvement of Trajectory Classification With auxiliary Information


Jérémie Leguay

Lorenzo Maggi

Moez Draief

Stefano Paris

Symeon Chouvardas

Algorithmes en ligne pour le contrôle d'admission

Keywords: Réseaux logiciel, Algorithmes en ligne, Apprentissage, Routage, Contrôle d'admission

En délégant le contrôle d?un réseau à des serveurs de calcul puissants, les architectures orientées logiciel permettent désormais l'utilisation en ligne d'outils apprentissage et d'optimisation. Dans ce contexte, ce papier étudie le problème du contrôle d'admission pour maximiser dans le temps le volume de trafic accepté. Il formule le problème de routage comme un problème de sac à dos, puis propose des adaptations d?algorithmes en ligne issues de la littérature. Il suggère ensuite d'utiliser l'algorithme expert intitulé Strategic Expert meta-Algorithm (SEA) pour décider en temps réel de l'algorithme en ligne à utiliser en fonction des conditions de trafic. Le papier présente une évaluation de cette architecture dans des conditions réalistes et illustre le gain en terme de trafic réseau accepté.

Introduction

Software-Defined Networking (SDN) technologies provide programmable data planes that can be configured from a remote controller. This separation between control and data planes creates an opportunity to implement routing processes that are more efficient than classic ones: in fact, the controller can take real-time decisions at a (logically) centralized location using an accurate and global view of the network.

A key task of the SDN controller is the Admission Control (AC) on incoming connection requests. Its goal is to gracefully manage service requests when the network becomes loaded. AC accepts or drops new requests depending on the resource availability. Non-myopic decisions have to be made with the aim of maximizing a given profit, such as the total accepted throughput, the financial revenue or the quality of service experienced by users. Nowadays, most of the deployed AC procedures are threshold-based. They use max-, min-, exclusive-and non-exclusive-limits on resource portions that the network operator can define for different classes of flows. The main problem here lies in defining the threshold in a dynamic fashion, as the optimal configuration depends on the network traffic conditions, which fluctuate over time.

In this paper we wish to raise the awareness that the ability of SDN controllers to centrally manage the network is an opportunity to revisit the way AC is performed. More specifically, we propose to extend online algorithms originally conceived for covering and packing problems. We then take a step further, and we propose to exploit the computational power offered by SDN controllers to implement machine learning techniques to boost performance (e.g., in terms of accepted throughput) via expert meta-algorithms, which are able to adaptively track the best AC algorithm without knowing the traffic statistics a priori. Specifically, we pinpoint a meta-algorithm called Strategic Expert meta-Algorithm (SEA) [START_REF] Pucci | How to combine expert (and novice) advice when actions impact the environment?[END_REF] which shows theoretical guarantees under our reactive scenario, on which there is little research.

Offline Admission Control

We represent the network as a capacitated graph G(V, E), where V and E are the set of nodes and directed edges in the graph, respectively. Let n = |V | denote the number of nodes in the graph. Each link e ∈ E has capacity u e . Connection requests arrive sequentially, and we denote the set of all connection requests by K.

The i-th request, with guaranteed bandwidth † r i , is described by a source-destination pair (s i , d i ), a pair of non-negative starting and ending times (t s i ,t f i ) and a profit b i . We denote P i as the set of feasible paths for connection request i (if i is accepted). We define f (i, p) as the portion of flow i that has been allocated to path p ∈ P i . Since we do not deal with fractional routing, f (i, p) ∈ {0, 1}, i.e., only one path can be used for each flow. The objective of the offline admission control problem is to maximize the total profit over the whole sequence of requests, known a priori, as follows:

max f ∑ i∈K ∑ p∈P i b i f (i, p) (1)
s.t. ∑ i∈K ∑ p∈P i |e∈p f (i, p) r i (t) ≤ u e , ∀ e ∈ E, t ≥ 0 (2) ∑ p∈P i f (i, p) ≤ 1, ∀ i ∈ K f (i, p) ∈ {0, 1}.
Nevertheless, solving (1) is not possible in practice: the controller receives information on the arrival and departures of requests as soon as they occur, and it has to make a decision on-the-fly.

Online Admission Control Algorithms

Traditionally, online algorithms for admission control fall into two categories: i) worst-case and ii) averagecase. i) Worst-case algorithms are characterized by max-min performance guarantees under specific worstcase scenarios where a malicious adversary chooses the worst possible sequence of connection requests.

Worst-case Admission Control (AC) Algorithms

Among the worst-case scenario AC algorithm, we first mention AAP algorithm [START_REF] Baruch Awerbuch | Throughput-competitive on-line routing[END_REF], taking admission control decisions based on the current utilization of network links. It computes path costs over a modified network graph where weights depend exponentially on the link utilization. This trick aims at pre-emptively driving traffic away from the links being highly utilized. The acceptance decision is based on a comparison between the cost of accepting the request and the resulting maximum future accepted throughput. Authors of AAP showed that the algorithm has a competitive ratio of O(log(n)) for any sequence of relatively small requests. This means that the number of accepted requests is in the worst-case O(log(n)) smaller than the number of requests that could be routed by the optimal solution of the offline problem in (1). AAP may not be easy to implement in reality, as it requires the a priori knowledge of requests duration. However, Buchbinder et al. proposed in [START_REF] Buchbinder | Online primal-dual algorithms for covering and packing[END_REF] a primal-dual framework to derive a practical algorithm with the same performance guarantees. The rationale behind it is that, when a demand arrives, the corresponding primal and dual variables are set while maintaining feasibility in both problems and while making sure that the derivative of the primal objective subject to the new dual variable evolves linearly with respect to primal variables, as proposed by [START_REF] Buchbinder | Online primal-dual algorithms for covering and packing[END_REF]. The second constraint guarantees the competitiveness of the algorithm. In the same manner as AAP, the acceptance decisions is taken by comparing the request cost (primal cost increase) and its profit (dual cost increase). We describe the steps of the Primal-Dual version of the AAP algorithm in Alg. 1, by using more efficient incremental updates of the primal variable x e .

Beyond Worst-case AC Algorithms

We now turn our attention towards online average-case (also called "stochastic") algorithms, showing good expected performance under random traffic conditions. Agrawal et al. [START_REF] Agrawal | Fast algorithms for online stochastic convex programming[END_REF] have proposed a fast algorithm with multiplicative updates to solve this issue. This recent algorithm described in Alg. 2 works for general convex problems. It applies to i.i.d. and random order inputs. It solves an online convex problem where the objective function is defined as the difference between the sum of rewards and the cost of accepted objects. The updates of θ and w are standard multiplicative weight updates used in the context of Algorithm 1 Primal-Dual AAP Algorithm [START_REF] Buchbinder | Online primal-dual algorithms for covering and packing[END_REF] Initialize x e = 0 function ROUTE(request j) if ∃ a path P ∈ P j of cost < 1 in the graph weighted by x e then Route request j on P for each edge e ∈ P do x e = x e exp ln(1+n).r j ue + 1 n (exp ln(1+n).r j ue -1) end for else Reject request j end if end function online optimization. Moreover, Agrawal et al. provide an alternative definition of competitive ratio as the ratio between the average reward and the average optimal reward. In this sense, the algorithm is 1 -O(ε)competitive for any ε > 0 such that min(B, k.OPT) > log(|E|)/ε 2 for a simple online covering packing problem. We adapted it to our multi-commodity flow problem.

The Primal-Beats-Dual (PBD) algorithm has also been introduced by Kesselheim et al. in [START_REF] Kesselheim | Primal beats dual on online packing lps in the random-order model[END_REF] for online packing problems with a finite number of objects. We considered it in our evaluation but it to solve a large LP at each step, which impacts on its scalability.

Admission Control with Experts in SDN

As shown in Sec.5, there is no algorithm that outperforms all the other ones under all traffic conditions. Hence, due to the unpredictable nature of traffic, it proves difficult to know a priori the identity of the best algorithm to be utilized. We thus need an online meta-algorithm that, based on past decisions and past rewards obtained by the different AC algorithms, can track and follow the best AC algorithm in hindsight.

This setting is classical in machine learning, and it is called prediction with expert advice. To be more formal, let us define m i, j as the traffic volume accepted by the AC algorithm i when request j arrives. The meta-algorithm takes its decisions iteratively, based on the profit m i, j obtained by each AC algorithms i over past decision instants j ≤ j. The bulk of the literature focuses on proving theoretical performance bounds in the basic non-reactive scenario where the action taken by the decision maker does not affect the state of the system. Nevertheless, our AC scenario is clearly a reactive one, since the decision taken at time t also influences the decisions (and the profits) of the AC algorithms at future time instants t > t.

To this aim, we hence propose to use Strategic Expert meta-Algorithm (SEA), described in Alg. 3. SEA selects algorithm i for an increasing number N i > 1 of consecutive steps, and does not revisit its choice at each new connection request. This allows each online AC algorithm to approach its asymptotically average performance. Moreover, SEA is only based on the profit effectively obtained by each algorithm when it has been actually selected. In other words, SEA does not exploit the information in hindsight on the profits that would have been obtained if a different algorithm had been used (as for FLA).

SEA's performance guarantees are evaluated in terms of the regret with respect to a (non implementable) oracle which steadily selects the algorithm with best average performance, known beforehand. In our

Algorithm 2

 2 Agrawal's Algorithm[START_REF] Agrawal | Fast algorithms for online stochastic convex programming[END_REF] Initializeθ 1,e = 1 1+|E| , ∀e ∈ E. Initialize w 0,e = 1, ∀e ∈ E. Initialize Z = OPT (B/T )function ROUTE(request j) Consider G(V,E) with edge cost of Zθ j,e , ∀ e ∈ E if p is a feasible min cost path in G then Route request j on p Perform the following multiplicative updates: w j+1,e = w j,e (1 + ε) (r j -u e /T ) , ∀ e ∈ E AND θ j+1,e = w j,e 1+∑ k w j,k , ∀e ∈ E else Reject request j end if end function

† For simplicity, we will denote r i (t) = r i for all t ∈ [t s i ,t f i ] and r i (t) = 0 otherwise.

Algorithm 3 Strategic Expert Meta-Algorithm (SEA) [START_REF] Pucci | How to combine expert (and novice) advice when actions impact the environment?[END_REF] Set M i = N i = 0 for each expert i. Set k = 1. function EXPERT SELECTION(Sequence of requests)

With probability 1/k perform an exploration phase, namely, choose an expert i from the uniform distribution over 1, . . . , N; otherwise, perform an exploitation phase, namely, choose an expert i with maximum M i . (If such i is not unique, select one out of a uniform distribution)

Set N i = N i + 1. Follow expert i for the next N i requests. Denote by R the average payoff accumulated during the current phase (i.e., these N i stages), and set

reactive scenario and under some stationary conditions on the system (in our case, of the traffic load on the links) SEA performs on average and asymptotically as well as the oracle ([dFM03], Thm. 5.1).

Performance Evaluation and Conclusion

To evaluate the algorithms under realistic conditions, we used a dataset captured in 2006 by Uhlig et al. on GEANT. We evaluated the percentage of rejected demands. The traffic matrix is generated with Poisson arrivals at rate λ demands/s. Demands have an equal size of 200 Mb/s and a duration exponentially distributed with mean 30s. We considered random source-destination pairs. We compared the results with the Greedy policy, that accepts all requests on the minimum cost path whenever there is enough capacity.

Remarkably, we notice that there is no algorithm that outperforms all the others under all traffic conditions, although their outperform Greedy This behavior naturally calls for a machine learning technique able to track the best AC algorithm under unknown traffic scenarios.

In Fig. 1 we show the performance of the SEA expert meta-algorithms which runs on top of three computationally efficient online algorithms: Greedy, AAP-pd and Agrawal's AC algorithm. As expected, SEA always outperforms the naive FLA.

The ultimate aim of this paper is to raise the awareness that the new centrally managed SDN architecture, combined with the computational power of the SDN controller, calls for a revisit of admission control algorithms, studied in the computer science literature but never really implemented in practice.


R Bertoni

M Lorenc

T Graber

R Henning

K Moffat

J.-F Létard

E Collet

Cooperative elastic switching vs laser heating in [Fe(phen) 2 (NCS) 2 ] spin-crossover crystals excited by a laser pulse

Spin-crossover crystals show multi-step responses to femtosecond light excitation. The local molecular photo-switching from low to high spin states occurs on sub-picosecond timescale. It is followed by additional conversion due to elastic (ns) and thermal (µs) effects. In [Fe(phen)2(NCS)2] crystals discussed herein, the thermal switching can be made unobtrusive for the investigation of cooperative elastic switching. We evidence a cooperative transformation induced by lattice expansion through elastic coupling between molecules in the crystal, where up to 3 molecules are transformed per photon.

Introduction

The behaviour of molecules in solid state differs from molecules in solution because of stronger interactions mediated through the lattice. One of the most striking examples is the existence of cooperativity where the responses of elementary components of the material add up in a non-linear fashion. This is especially true for photoinduced phase transitions, where the excitation by a light pulse can abruptly generate cooperative transformations towards a new macroscopic state. [START_REF] Nasu | Photoinduced Phase Transitions[END_REF] Such transformations can be triggered in an ultrafast way where a femtosecond laser pulse impacts the material. During photoinduced phenomena, several microscopic and macroscopic degrees of freedom play their part on significantly different time and length scales. For example, transformations at the molecular scale can be coherently activated [START_REF] Ishikawa | [END_REF]3,4 on the timescale of hundreds of femtosecond. Slower long-range structural reorganizations may induce conductivity or ferroelectricity for example. 5,6 In addition to localized processes, the entire material can react through lattice deformations and volume expansion occurring on the acoustic time scale 7,8,9,10

. The understanding of the response of the material at this time scale is of paramount importance because it may stabilize the newly formed photoinduced phase. Here we focus our attention on Spin-crossover (SCO) molecular materials, which are prototypical photoactive materials showing molecular bistability between electronic states differing in physical parameters such as volume, colour, magnetic susceptibility, etc... [START_REF] Halcrow | Spin-crossover materials[END_REF] For an Fe(II) system like the [Fe(phen) 2 (NCS) 2 ] single crystal investigated here, the low spin (LS) S=0 state corresponds to a t 2g . In the less bonding HS state the <Fe-N> bonding is elongated by 10% and the molecule swells. The SCO systems are known to be photoswitchable between LS and HS states. [START_REF] Halcrow | Spin-crossover materials[END_REF] When induced in the solid state by a femtosecond laser pulse, this process is the initial trigger of a complex out-of-equilibrium dynamics, where different degrees of freedom act at different time and length scales. The qualifying term "molecular materials" unfolds then its genuine meaning: the material response to light stimuli far exceeds the sum of the individual responses of constituting molecules. We have recently demonstrated the possibility of taking advantage of this effect to generate through inter-molecular elastic coupling an elastically-driven cooperative response. Our previous studies have shown that the out-of-equilibrium transformation of SCO crystals triggered by a fs laser pulse involves three main steps (Fig. 1). [13][14][15][16] The light pulse locally photoswitches molecules from LS to HS via intersystem crossing through metal-toligand charge-transfer excitation (MLCT, t 2g 5 e g 0 L 1

). It was shown by ultrafast techniques that during this photoswitching step, the molecular breathing that accompanies the LS-to-HS photoconversion occurs in the solid within less than 200 fs whereas vibrational cooling occurs on the ps timescale. 4,[17][18][19][20][21] An initial fraction of HS molecules X HS h is thus photo-generated in the crystal. A second increase of the HS fractions up to X HS El occurs later on the ns timescale. The internal lattice pressure due to HS molecular swelling and the lattice heating induces lattice expansion, driving this additional switching process: the elastic coupling between molecules favours the HS state of higher molecular volume. Because of the global temperature increase of the crystal (10s K) a third "thermal" step sets in when the HS population thermally equilibrates to X HS Th with the surrounding hot lattice, typically within 10 µs. The HS molecular states generated by pulsed laser excitation are transient when induced from the low temperature LS state (typically in the ms range above 100 K). At very low temperature the lifetime of photoinduced states is much longer and complete conversion can be easily reached by weak cw laser excitation. [START_REF] Halcrow | Spin-crossover materials[END_REF] Also, a single laser pulse can drive a complete conversion of the material inside thermal hysteresis.

22,23

As both elastic or thermal effects play their part in the out-ofequilibrium process, it is important to separate these two contributions. Here we investigate the response of photoexcited [Fe(phen) 2 (NCS) 2 ] single crystals by combining time-resolved optical spectroscopy, probing the change of electronic state (HS fraction), and time-resolved x-ray diffraction, probing crystalline lattice change. The first-order nature of the phase transition at thermal equilibrium allows disentangling elastic and thermal effects.

Results and discussion

Thermal transition and the photo-induced steady state

The thermal phase transition of [Fe(Phen) 2 (NCS) 2 ] crystals was intensively investigated by different spectroscopic and scattering techniques. [24][25][26][27][28] In this study we used single crystals of typical size 200µm×200µm×30 µm 3 , which undergo a first-order phase transition from purely LS phase below 180 K to purely HS phase above. The discontinuous change from LS to HS states is accompanied by a large structural reorganization both at the molecular level around the Fe-N 6 core and the macroscopic scale. Fig. 2 shows data reported by Cammarata et al, 4 which agrees with previous crystallographic studies. 26,[START_REF] Gallois | Structural Changes Associated with the Spin Transition in Fe(phen)2( NCS)2: A Single-Crystal X-ray Investigation[END_REF][START_REF] Granier | High-Pressure Single-Crystal X-ray Diffraction Study of Two Spin-Crossover Iron( II) Complexes: Fe(Phen)2(NCS)2 and Fe(Btz)2(NCS)2[END_REF] The average bond length <Fe-N> elongation from 1.97 Å in the LS state to 2.16 Å in the HS state, is characteristic of Fe(II) SCO materials. Through the LIESST effect (Light Induced Electronic Spin State Trapping) discovered in the 80's by Hauser and co-workers, 31 it is also possible to reach long-lived HS state by weak cw light excitation at low temperature. Photo-crystallography studies performed at 25 K indicate a global conversion of the crystal after excitation at 660 nm with a cw laser, associated with an elongation of <Fe-N> up to 2.16 Å and confirming previous studies on the same compound after excitation at 647 nm. 26 An important feature of the LIESST phenomenon is that its quantum efficiency is close to 1, i.e. almost every absorbed photon switches a molecule from LS to HS state. From the macroscopic point of view, the increase of molecular volume between LS and HS states induces an important lattice expansion, both during the thermal and the photoinduced LS to HS states conversions (Fig. 2). 4,9,[START_REF] Bertoni | [END_REF]17 This can be seen in the discontinuous evolution of the lattice unit cell volume. Figure 3 shows the evolution of the intermolecular distance due to molecular swelling between LS and HS states. In addition to structural changes, another consequence of the change of electronic state during the spin state transition is the change of optical properties. The HS and LS states have different absorption bands due to electronic and structural reorganizations. [4][5][6][7][8][9][10][START_REF] Halcrow | Spin-crossover materials[END_REF][START_REF] Bertoni | [END_REF][13][14][15][16][17][18][19][20][21] This provides a spectroscopic probe of the spin state switching, often used to characterise SCO phenomena. [START_REF] Halcrow | Spin-crossover materials[END_REF] Here we monitor the thermal transition in [Fe(Phen) 2 (NCS) 2 ] single crystals through optical transmission measurements (Fig. 2b) at 950 nm. In this spectral zone, the optical transmission decreases in the HS state of higher optical density than the LS one. The photoexcitation of purely LS crystals by a femtosecond laser pulse into the MLCT band (in the 640-660 nm range) efficiently switches molecules from LS to HS states. 4,21 The time-resolved optical experiments developed at the Institut de Physique de Rennes allow monitoring the optical transmission change after femtosecond photoexcitation, continually from 100 fs to ms, and thereby covering the entire timescale of non-equilibrium dynamics. The time-resolved optical transmission change at a probing wavelength of 950 nm is then scaled in terms of HS fraction X HS . In addition, for investigating the role of lattice expansion in the process, we also performed timeresolved x-ray diffraction studies at the BioCARS beamline 32 at the Advanced Photon Source, Argonne National Laboratory (see experimental section).

Response of a single crystal impacted by a femtosecond laser flash

We photoexcited 25 µm thick [Fe(Phen) 2 (NCS) 2 ] single crystals in the tail of the LS MLCT band with a 650 nm light pulse. The penetration depth at this wavelength is on the order of the thickness of a single crystal. The evolution of the HS fraction X HS after fs laser excitation in the LS phase at 140 K, measured by optical spectroscopy, is shown in Fig. 4 for different excitation densities. The experiment is conducted 40 K below the phase transition to avoid any residual thermal effect and to allow the complete recovery of the LS phase within 1 ms. Our experimental data clearly show a sequence of three steps, similar to the ones already reported in Fe(III) spin-crossover crystals and summarized in Figure 1. [START_REF] Bertoni | [END_REF][13][14][15][16] First, the absorption of light at the molecular level locally photo-switches a small fraction X HS h of molecules (typically 2% here) from LS to HS states. We have thoroughly investigated this process in [Fe(Phen) 2 (NCS) 2 ] crystals with femtosecond optical and x-ray spectroscopies 4,21 and established that it occurs via ultrafast intersystem crossing. This process involves local coherent structural dynamics and is accompanied by a molecular swelling within 200 fs. X-ray diffraction experiments (Figure 5), performed at 140 K for an excitation density of 20 µJ/mm 2 , reveal that during this initial photoswitching process at the molecular level, the lattice has no time to expand. Indeed the values of the unit cell parameters remain constant up to the ns time scale. The lattice expansion is observed only after a few ns, and occurs concomitantly with a second increase of the HS fraction, reaching X HS El after 20 ns (Figure 4). The microscopic crystal expansion results from the establishment of a mechanical equilibrium within the crystal: the photoswitched molecules of higher volume exert a negative (or internal) pressure on the lattice, thus causing its expansion. Such effect stems from elastic properites and is responsible for the selfamplified molecular transformation, as recently demonstrated both in experiment and theory. [START_REF] Bertoni | [END_REF] It should be underlined that the cell volume expansion is limited by the propagation of strain waves through the sample. Given the typical acoustic wave velocity (~2000 m.s -1

) in these materials and the anisotropic shape of the crystals used (lozenge plate) the expansion occurs faster on the short crystal dimension (25 µm along the c axis) than on the long dimensions (200-300 µm along the a,b axes). At excitation densities above 40 µJ/mm 2 there appears a third step in which the HS fraction reaches X HS Th on the 10 µs regime (Figure 4). This effect is associated with the thermal population of the HS state, due to equilibration of the SCO molecules with the hot crystalline lattice, as the absorbed photon energy is converted to heat. As the HS state is more and more populated, the lattice expansion increases anew. The response on the photo, elastic and thermal steps to different excitation densities (expressed in number of incident photon per 100 molecules) is shown in Fig. 4b. The photoswitched fraction X HS h depends linearly on the laser fluence.

By comparing the number of incident photons per molecules in the pump laser pulse with the fractions of photo-switched molecules, we estimate that essentialy every photon switches one molecule. This confirms the quantum efficiency close to unity observed in SCO solids. 4,[START_REF] Bertoni | [END_REF] The photo-response on the elastic step is between 2.5 and 3 times larger than on the photoswitching step. A linear fit of the fraction of HS molecules on the elastic step with the number of incident photons indicates that in average 1 photon transforms 2.66 molecules (green line in Fig. 4b).

The photoresponse on the thermal step shows a non-linear response and strongly increases with excitation density. To unmbiguously assign the last step to thermal effects, we performed similar optical spectroscopy experiments at 140 K and 160 K with an excitation density fixed at 30 µJ/mm 2 (Figure 6). At 140 K only the photoswitching and elastic steps are observed resulting in a twostep response in X HS . At 160 K the last thermal step is now observed and characterized by a higher HS population in the µs domain. It is accompanied by a weak lattice expansion in the 100 µs regime observed at 160 K under the same 30 µJ/mm 2 excitation density (orange circle, Fig. 5). The 2 Å 3 volume expansion after the thermal step, with respect to the equilibrium values at negative delays, can be compared to the thermal expansion observed at thermal equilibrium. As shown in figure 5, before laser excitation the volume increase between 140 K and 160 K is

2 Å 3 (from 2206 to 2208 Å 3
). It is therefore reasonable to estimate from the 2 Å 3 volume expansion at 100 µs delay, that 30 µJ/mm 2 excitation density generates a 20 K laser heating. This means that in the time domain the lattice temperature approaches 180 K i.e. the thermal phase transition temperature. It explains the significant thermal population observed at 160 K and for 30 µJ/mm 2 . In contrast photoexcitation at 140 K with the same (Fig. 6) or lower (Fig. 4 &5) laser fluence is not enough to populate the HS state thermally because of the first-order nature of the phase transition, as the HS fraction changes discontinuously from 0 to 1 around 180 K.

In Figure 6, which compares the photo-reponse after identical laser excitations at 160 and 140 K, it appears that up to 1 µs the two curves are identical within the experimental accuracy. At 140 K, the HS molecules converted after the elastic step return to the LS ground state within 1 ms and no further increase of the HS fraction is observed. At 160 K, a significant thermal population occurs on the 10-100 µs scale. The difference (top panel, black triangle, Fig. 6) between the photoexcitation at 160 K and 140 K reveals only the thermal population of HS and demonstrates that this process occurs in the µs-ms regime. The experiment performed at 140 K for different laser fluence also confirms the thermal nature of the µs conversion (Fig. 4). For excitation density below 30 µJ/mm 2 there is no thermal effect because the temperature jump does not bring the system close enough to the phase transition to populate the HS state. Higher excitation densities induce a larger heating effect. For example at 60 µJ/mm 2 an almost 40 K temperature jump occurs, if we assume a linear dependence of the temperature jump with the deposited energy on the system. This is enough to bring the system close to the phase transition and populate the HS state thermally. This result is different from previous studies on weakly cooperative Fe(III) SCO materials, 9,10,17 for which thermal conversion is detected because unlike in a first-order transition, the crossover spans noncritically a broad temperature range. In the present case, it is possible to avoid thermal population when the experiment is performed at a low temperature, or with a low excitation density.

Conclusions

In the out-of-equilibrium process following femtosecond laser excitation of SCO single crystals, three main steps were identified, with characteristic timescales: photo-switching (ps), elastic switching (ns) and thermal switching (µs). One extraordinary aspect of such multiscale response of the material is that it involves different physical processes: from quantum physics and intersystem crossing during the photo-switching step, mechanics during the elastic step and thermodynamics during the thermal step. Our study of [Fe(II)(Phen) 2 (NCS) 2 ] SCO material corroborates the previous findings and shows the universality of this picture in SCO materials. The combination of time-resolved x-ray diffraction and optical spectroscopy reveals clearly the thermal nature of the last step, which can be avoided with appropriate excitation density or initial temperature. In this cooperative SCO material we observe a self-amplified response during the elastic step after excitation of a purely LS state where up to 3 molecules are transformed per photon. Compared to isolated photoactive molecules, with a quantum efficiency attaining at best unity, these results show that mechanical forces induced by absorption of light can induce cooperative switching in volume-changing materials and surpass the single molecule limits. Our previous studies on Fe(III) systems 12 evidenced that the elastic response appears only above a threshold excitation density. It is not the case here as the elastic response varies linearly with excitation density (Fig. 4b). This difference may be attributed to the larger molecular deformation in the Fe(II) systems, where the Fe-ligand distance elongates by 0.2 Å, whereas for Fe(III) it does only by 0.15 Å. The energy cost or the elastic coupling is therefore larger as underlined by the first-order character of the phase transition of this compound, compared to the more gradual conversion in Fe(III) systems studied previously. This cooperativity may favor local clustering around photo-switched molecules, resulting in a more linear response. This point will be investigated in future with the mechano-elastic model applied to rationalize this complex out-of-equilibrium transformation at the material scale.

12

Experimental methods

Time-resolved x-ray diffraction and optical transmission measurements were performed on single crystals, with typical dimensions of (250±50)(200±50)(30±5) µm

3

. The laser was propagating along the 30 µm thickness corresponding to the crystalline axis c. The single crystals were mounted in nitrogen-flow cryostreams. The optical pump-probe studies were performed at the Institut de Physique de Rennes with the experimental setup described in detail in ref. 9. Single crystals were excited with a 80 fs laser flash and their optical transmission was probed by a delayed 80 fs laser probe. The delay between pump and probe was controlled by combining a mechanical translation stage to adjust the optical path difference for sub-ns measurements with an electronic synchronization of the pump and probe amplifiers for measurements in the ns to ms time domain. The experiments were performed with 500 Hz laser pump repetition rate, set to 650 nm on the MLCT band where it efficiently induces LS-to-HS conversion of the[Fe(phen) 2 (NCS) 2 ] crystal. The probe was set to 950 nm to monitor the LS-to-HS photoswitching dynamics through the variation of optical transmission, and from which X HS was extracted (see also ref. 14 and 21). Time-resolved x-ray diffraction experiments were performed at the BioCARS beamline at the APS synchrotron. 32 Single crystals were excited with a ps laser flash centered around 650 nm. Single x-ray pulses (15 keV) were selected by a fast chopper for probing the crystal at different delays, with a 40 Hz repetition rate. The 70 ps time resolution is limited by the xray pulses duration. Partial data were collected in transmission geometry at different pump-probe delays using a single sample rotation axis and the unit cell parameters were obtained from diffraction images with CrysAlis sofware. 33 The crystallographic data presented in Fig. 2 & 3 were obtained at the Institut de Physique de Rennes (see also ref. 4) on a four-circle Oxford Diffraction Xcalibur 3 diffractometer (MoK α radiation). The single crystals were mounted in an Oxford Cryosystems nitrogen-flow cryostat for experiments above 78 K and in an Oxford Diffraction Helijet cryostat for measurements at 25 K, where a 660 nm laser diode was used for photoexcitation.

Figure 1 .

 1 Figure1. Schematic drawing of the out-of-equilibrium dynamics in SCO crystals. The laser pulse locally photo-switches molecules from LS (blue) to HS (red) states within less than 1 ps (step 1). The molecular swelling and lattice heating induce lattice expansion driving cooperative elastic switching in the ns time window (step 2). The HS state is then thermally populated (µs) and accompanied by another crystal expansion (step 3).

Figure 2 .

 2 Figure 2. Fingerprints of the conversion between LS and HS states. a) Temperature dependence of the lattice parameter c, the volume V and the <Fe-N> bond length. At 25 K the photoinduced HS state is reached by excitation around 650-660 nm. b) Optical transmission at 950 nm.

Figure 3 .

 3 Figure 3. The molecular swelling between LS and HS states induces lattice expansion, shown in along the a crystalline axis.

Figure 4 .

 4 Figure 4. a) Response to femtosecond laser excitation of LS [Fe(Phen)2(NCS)2] at 140 K for different excitation densities "h", "El" and "Th" denote respectively the photoinduced, elastic and thermal steps. b) Evolution of XHS h (measured at 10 ps) scaling with the number of incident photon per 100 molecules. On the photo-switching step 1 molecule is transformed per photon, whereas on the elastic step (XHS El ) 2.5-3 molecules are transformed per photon. Heating effect appears at high excitation density (close to sample damage).

Figure 5 .

 5 Figure 5. a) Lattice response to femtosecond laser excitation of LS [Fe(Phen)2(NCS)2] at 140 K for an excitation density of 20 µJ/mm 2 and at 160 K for an excitation density of 30 µJ/mm 2 .

Figure 6 .

 6 Figure 6. Evolution of the HS fraction XHS after femtosecond laser excitation of LS [Fe(Phen)2(NCS)2] at 140 and 160 K for 30 µJ/mm 2 excitation density.

Acknowledgements

This work was supported by the CNRS, the Institut Universitaire de France, Rennes Métropole, the ANR (ANR-13-BS04-0002) and Europe (FEDER). K.M., R.H., T.G. and E.C. thank the University of Chicago for funding through the "France and Chicago Collaborating in the Sciences" (FACCTS) program. The BioCARS facility is supported in part by the Department of Energy (DE-AC02-06CH11357) and in part by the National Institutes of Health grant GM111072 to K.M.. The timeresolved set-up at BioCARS was funded in part through a collaboration with Philip Anfinrud (NIH/NIDDK).


Mahmoud Soua
email: mahmoud.soua@esiee.fr

Rostom Kachouri
email: rostom.kachouri@esiee.fr

Mohamed Akil
email: mohamed.akil@esiee.fr

Efficient Multiscale and Multifont Optical Character Recognition System based on robust Feature Description

Keywords: OCR System, Multiscale, Multifont, Feature Extraction, Feature Matching, SAD technique

Optical Character Recognition (OCR) is the process of translating images of text into a comprehensible machine format. Generally, an OCR system is composed of binarization, segmentation and recognition stages. Given an extracted binary character, the recognition stage ensures its description and decides its corresponding ASCII code. In this paper, we propose a new OCR system that aims to high speed, Multiscale and Multifont character recognition. Our proposal is based essentially on robust description using a new Unified Character Descriptor (UCD). In addition, a character type-face and fontsize recognition is performed to choose the adequate template for faster matching process. Obtained OCR Accuracy of our proposed System is 1.5x higher then that reached by Tesseract on the LRDE dataset.

I. INTRODUCTION

Optical Character Recognition (OCR) deals with the problem of recognizing optically processed characters. OCR systems involve three major stages to completely recognize characters; Binarization, Segmentation, and Recognition stages [START_REF] Prasad | Template matching algorithm for Gujrati Character Recognition[END_REF]. Firstly, Binarization separates the text characters from the background [START_REF] Soua | A new hybrid binarization method based on Kmeans[END_REF]. Then, Segmentation stage aims to locate text regions in the processed documents [START_REF] Kopf | Robust Character Recognition in Low-Resolution Images and videos[END_REF]. Finally, the Recognition stage consists on a very sensitive character description and decision steps.

Generally, state-of-the-art methods in the description step are based either on Template Description (TD) or on Feature Description (FD) [START_REF] Prasad | Template matching algorithm for Gujrati Character Recognition[END_REF]. In TD Methods [START_REF] Prasad | Template matching algorithm for Gujrati Character Recognition[END_REF][2], characters are described based on their pixel information. Despite their simple use, these methods are not robust on noisy characters. In addition, calculation are burden with unnecessary pixel description. Less complex, the FD methods [START_REF] Aghdasi | Automatic Licence Plate recognition system[END_REF][4] [START_REF] Farhad | An efficient Optical Character Recognition Algorithm using Aritificial Neural Network by Curvative Properties[END_REF] perform the description of characters based on some specific Features.

In the other hand, the decision stage is made according to a matrix or extracted feature matching. As mentioned above, pixel comparison in matrix matching is extremely sensitive to noisy characters [1][2]. This stage can be performed using classification methods such as RNA, SVM [START_REF] Aghdasi | Automatic Licence Plate recognition system[END_REF][4] [START_REF] Farhad | An efficient Optical Character Recognition Algorithm using Aritificial Neural Network by Curvative Properties[END_REF]. They give interesting results, however they still complex comparing to matching process. Indeed, simple feature matching ensures a good trade off between accuracy and computation complexity [START_REF] Farhad | An efficient Optical Character Recognition Algorithm using Aritificial Neural Network by Curvative Properties[END_REF].

In this paper, we propose a novel Multiscale and Multifont OCR System based on a robust feature description. Firstly, we ensure the right template selection. Then, we compute a Unified Character Descriptor (UCD) and a fast Matching process is performed.

In the following, we present our proposed OCR System in section 2 where the new template selection and Character Recognition based on UCD (CR-UCD) method is explained. The obtained results are shown and discussed in section 3. We design a complete OCR System that performs character recognition on input image containing text information. As seen in Figure 1, our OCR System includes the stages of binarization, segmentation, template selection and recognition. In the binarization stage we use our recently proposed Hybrid Binariaztion Based on Kmeans (HBK) [START_REF] Soua | A new hybrid binarization method based on Kmeans[END_REF][13] method to separate correctly the text and the background even on noisy images. Then, binarized image is segmented with the Page Layout Analysis (PLA) [START_REF] Smith | Hybrid Page Layout Analysis via Tab-Stop Detection[END_REF] part of the well known Tesseract 3.02 engine [START_REF] Smith | An Overview of the Tesseract OCR Engine[END_REF]. In the template selection stage, character font-size and type-face are recognized and subsequently the adequate template is chosen. Finally, we perform our CR-UCD recognition in which each character is represented firstly in the Description stage using one UCD feature. Then, characters are recognized based on feature matching with the selected template using the Sum of Absolute Difference (SAD). An example of letter 'A' recognition with our proposed method is given in Figure 2, we assume that the considered letter 'A' is extracted from the word 'Académies'. In this figure, the different proposed stages of our new method, namely: the Template Selection, Description and Decision stage are illustrated. In the first stage a Template Selection Descriptor (TSD) is generated to allow the appropriate template choice. In an other hand, the description stage extracts the Unified Character Descriptor (UCD) from each character and feeds it into the decision stage. In this stage the extracted UCD vector is matched with the Fig. 2: Letter 'A' recognition with the proposed method closest character in the selected template. These stages are given with more details in subsequent sections. Firstly, we explain the Description stage. Then, we discuss the template selection and we present the decision ones.

A. The Description Stage

In this section, we present our Feature Extraction Strategy. The aim of this stage is to employ a sufficient number of characteristics that helps to discriminate characters efficiently. For this, we perform in our work a Character Segment Extraction by using a simple Edge Detection. We extract then horizontal and vertical character segments. In the Horizontal Edge (HE) extraction, for each pixel P ij we compute the right edge according to Equation [START_REF] Prasad | Template matching algorithm for Gujrati Character Recognition[END_REF]. Where i ∈ [1..h] and j ∈ [1..w], given that h and w are respectively the height and width of the character bounding box. We note by 0 a black pixel and by 1 a white one.

P i,j ∈    HE If (P i,j = 0) And (P i,j+1 = 1)
HE If ((P i,j = 0 And (P i,j+1 = 0) Or (P i,j = 1))

(1) In the other side, the Vertical Edge (VE) extraction, for each pixel P ij is computed as shown in Equation [START_REF] You-Qing | A Recognition Method of Car License Plate Characters Based on Template Matching Using Modified Hausdorff Distance[END_REF].

P i,j ∈    VE If (P i,j = 0) And (P i+1,j = 1)
VE If ((P i,j = 0 And (P i+1,j = 0) Or (P i,j = 1))

(2)

Generally, when we re-scale characters or when we handle the same character font in many sizes, we change the character morphology and some distortions can appear in the character description. This issue is named the aliasing behavior [START_REF] Smith | Hybrid Page Layout Analysis via Tab-Stop Detection[END_REF]. As shown in Figure 3, due to this phenomena, one character can have a different number of segments on multiple Font-Sizes. For this, we propose to check the extracted segments to merge some edges based on their neighborhood. To ensure an efficient and scale invariant method, we use each time a variable neighbour pixel number for the check. This parameter is adjusted according to h and w values of the processed character (its size). This proposed technique allows to improve the description of multi-scale character.

As shown in Figure 2, the Character Segment Description (CSD) of each Horizontal and Vertical extracted segments is firstly ensured by the Character Segment Position (TMB and LMR). Then one Horizontal Merged CSD is constructed while adding respectively H-TMB and H-LMR of each Horizontal segment. To which we concatenate other features like the total Horizontal Character Segment Number (H-NBS), the Horizontal Character Barycentre Coordinates (H-Bx and H-By), and the Horizontal Character Ratio (H-R). We do the same thing for the Vertical Merged CSD computation except the Character Ratio (R) feature which is not used in this case. The employed features in this work are presented and explained with more details in the following subsections:

1) Character segment Number Feature (NBS): We consider that this feature increases the discrimination between characters that do not have the same number of segments.

2) Character Segment Position Feature (TMB, LMR): As shown in Figure 2, the CSD vector is composed of two features: The first one, TMB, refers to the Top, Middle and Bottom segment positions in the character bounding box. We compute TMB according to Equation (3).

T M B =        1 If (S y > h 2 ) (Top) 2 If (S y = h 2 ) (Middle) (3) 3 If (S y < h 2 ) (Bottom
LM R =        1 If (S x < w 2 ) (Left) 2 If (S x = w 2 ) (Middle) (4)
3 If (S x > w 2 ) (Right) with S x is the starting x coordinates of the profiled segment and w is the width of the character bounding box.

3) Character Barycentre Coordinate Feature (Bx, By): Starting from the motivation that different character shapes have different barycentre positions, as shown in Figure 4, we propose a novel and simple barycentre computation technique. Indeed, we consider the polygon composed only with the starting segment pixels. Then, we compute the x and y coordinates of the corresponding barycentre for both Horizontal and vertical Merged CSD Vectors.

Figure 5 shows a Multiscale representation of characters 'A' and 'J'. As we can see, the Ratio feature (R) does not changes for the same character with different sizes (Equation ( 6)):

   R( A ) = w00 h00 = w01 h01 = w02 h02 = w03 h03 = w04 h04 = w05 h05 R( J ) = w10 h10 = w11 h11 = w12 h12 = w13 h13 = w14 h14 = w15 h15 (6)
However, in the same scale, the Ratio feature (R) allows to discriminate easily the two different characters 'A' and 'J' (Equation ( 7)):

           Scale 0 : [R( A ) = w00 h00 ] = [R( J ) = w10 h10 ] Scale 1 : [R( A ) = w01 h01 ] = [R( J ) = w11 h11 ] . . . Scale 5 : [R( A ) = w05 h05 ] = [R( J ) = w15 h15 ] (7)
In fact, the Ratio feature (R) ensures a Multiscale invariance and increases the discrimination between characters in the

B. The Template Selection Stage

In our proposed method, the template selection stage is performed through type-face and font-size recognition. Indeed, even if the above UCD is scale invariant, it stills more powerful on small size ranges. To improve the Multiscale characteristic of our method, we propose to compute effectively the character font-size to choose the appropriate template to use. In addition, for an identified type-face character it is possible to make the OCR system handle a document with less effort. Following we detail our template selection proposal.

1) Font-Size Recognition (x-h): Text line documents are composed of three typographical zones: the ascender, the xheight and the descender zones [START_REF] Abdelwahab | Optical Font Recognition Using Typographical Features[END_REF]. As shown in Figure 7, these zones are delimited by four virtual horizontal lines, Ascender, x-height, Base and Descender lines. The x-height Fig. 7: Typography of Latin Alphabets zone is the height of the body of lowercase letters referring to the distance between the baseline and the x-height one of lower-case letters in a typeface. For Latin script the Font-Size recognition is relatively an easy task [18][14]. In our proposal, we identify the character size with the help of the computation of the connected components x-h. Hence, we apply the following formula shown in Equation ( 8) [START_REF] Souza | Automatic filter selection using image quality assessment[END_REF].

x-h = Baselinex-height line ( 8)

2) Type-Face Recognition (σ): Type-face recognition can give details on the structural and the typographical design of characters. Type-faces can be divided into two main categories: serif and sans serif. We mean by serifs the small features at the end of strokes within letters. Printed type-faces without serifs are known as sans serif. 'A' letter with serif and sans serif type-faces is shown in Figure 8. 9)).

σ = 1 N -1 N i=1 (x i -µ) 2 ; µ = 1 N N i=1 (x i ) ( 9
)
where x is the measured thickness in one character block, N is the number of character strokes and µ is the average of the N measured thicknesses; i ∈ [1, N ]. To do this, we divide the character image generated by the layout analysis stage into 3x3 blocs. As shown in Figure 9, we determine in each block the thickness as the median of non overlapping lines with the block edges. During this process, vertical (stem, hairline, etc) and horizontal strokes (cross stroke, cross bar, etc) are checked. Very long or short lines (such as brackets, apex) that do not give a significant indication regarding the thickness are eliminated thanks to the 3x3 bloc division. As illustrated in Figure 10, the measured stroke thickness Standard Deviation (σ) is approximately constant between different Font-Sizes however it varies considerably between different type-faces.

C. The Decision Stage

The final stage of our work consists on matching the UCD vector of the input character to the different UCD vectors of the selected template (see Figure 2). The matching is performed by using the Sum of Absolute Differences (SAD) technique [START_REF] Lengwehasarit | Probabilistic partial-distance fast matching algorithms for motion estimation[END_REF] distance process.

III. EXPERIMENTAL RESULTS

In this section, we evaluate firstly our proposed UCD descriptor regarding the robustness against the 'i' algorithm [START_REF] Sharstry | i" A novel algorithm for optical character recognition[END_REF] descriptor when dealing with Multifont and Multisize constraints. Then, we evaluate our CR-UCD system against the well known Tesseract 3.02 OCR engine [START_REF] Smith | An Overview of the Tesseract OCR Engine[END_REF] on the same previously mentioned constraints.

A. The UCD Descriptor Evaluation based on OCR Accuracy

To evaluate the OCR-accuracy of our proposed UCD descriptor we use Multiscale and Multifont computer generated alphanumeric images containing number, upper-case, lowercase and special characters. For comparison reasons, we assess also the recently proposed 'i' Algorithm [START_REF] Sharstry | i" A novel algorithm for optical character recognition[END_REF]. Table 1 shows a comparison between UCD and 'i' Algorithm descriptors using OCR Accuracy. In this evaluation, each character is matched with the same character font and size in the template. Obtained results demonstrate that our method outperforms the 'i' algorithm reaching an average of 99% of OCR accuracy on all assessed font-sizes and Type-faces. The main issues encountered by our description is the miss-recognition of similar characters like 'I' and 'l' in the different scales of the sans serif Type-face. In future work, lexical methods can be used in that case to enhance the performance of our description.

In Table 2 we illustrate the performance of UCD and the 'i' Algorithm descriptors using one single template defined by the serif type-face and the 90pt font-size for all compared character fonts and sizes. We can see that in this case the Small characters refer to the core paragraph text, medium characters refer to subtitles and large size characters are considered as the titles. Empirically fixed, the employed thresholds to distinguish font-sizes are given in Table 3. To show a relevant evaluation, we compare our proposed CR-UCD OCR system to the Tesseract [START_REF] Smith | An Overview of the Tesseract OCR Engine[END_REF] one. To make a fair comparison, we disabled the Tesseract dictionary option that we did not handle yet in our work. Figure 13 shows that the proposed CR-UCD method outperforms Tesseract one by reaching 1.5x higher OCR Accuracy on the 125 magazine documents of the LRDE Dataset.

IV. CONCLUSION AND FUTURE WORK

In this paper, we proposed an efficient OCR system based on a new and robust description. Thanks to the used UCD descriptor, the template selection, and the employed fast feature matching our proposal offers a high OCR Accuracy compared to well known description and recognition methods. In future work, we intend to improve the proposed OCR efficiency by using additional semantic rules to be able to differentiate similar characters like 'I' and 'l'. In addition, the template selection process could be enhanced to distinguish thick and italic characters. We plan also to use dictionary to enhance the final recognition results.

Fig. 1 :

 1 Fig. 1: Block diagram of our proposed OCR System

Fig. 3 :

 3 Fig. 3: Dependency of scale and character segment discontinuity: a. Input Images, b. Extracted Horizontal Edges in one scale, c. Extracted Horizontal Edges in another scale, d and e. A corresponding zoom views of the extracted segments

  ) with S y is the starting y coordinates of the profiled segment and h is the height of the character bounding box. The second feature is LMR. It refers to the Left, Middle and Right segment positions in the character bounding box. It is computed as shown in Equation (4).

4 )Fig. 4 :

 44 Fig. 4: Barycentre (B) computation process on different segment number letters: a. Starting pixels, b. Barycentre of polygon computation, c. Barycentre result

Fig. 5 :

 5 Fig. 5: Scale invariance of the character ratio feature of 'A' and 'J' letters with six different scales

Fig. 6 :

 6 Fig. 6: The Unified Character Descriptor (UCD): a. Character Segment Number Feature (NBS), b. Character Segment Position Feature (TMB, LMR), c. Character Barycentre Feature (Bx, By), d. Character Ratio Feature (R)

Fig. 8 :

 8 Fig. 8: 'A' letter with: a. Sans Serif and b. Serif type-faces

Fig. 9 :

 9 Fig. 9: Font recognition based on the splitting of characters in equal size blocks. Each block number refers to the median thickness of the considered non overlapping strokes

Fig. 10 :

 10 Fig. 10: Standard Deviation (σ) of stroke thickness in Multiscale and Multifont characters. a. Sans Serif with no/moderate transition in the strokes, b. Serif with dramatic thick/thin transition in the strokes

Fig. 11 :Fig. 12 :

 1112 Fig. 11: Sample of the LRDE-DBD documents

Fig. 13 :

 13 Fig. 13: OCR Accuracy evaluation of Tesseract and the proposed CR-UCD method

Table 1 :

 1 OCR Accuracy of UCD and 'i' Algorithm descriptors using Multifont and Multisize template

	Method	Type-face			Font-Size
			18	24	36	72	90
	'i' Algorithm	Serif Sans Serif	90 80	78 77	90 82	95 93	98 95
	CR-UCD	Serif Sans Serif	100 98	100 98	100 98	100 98	100 98

Table 2 :

 2 OCR Accuracy of UCD and 'i' Algorithm descriptors using one single font and size template

	Method	Type-face			Font-Size
			18	24	36	72	90
	'i' Algorithm	Serif Sans Serif	12 15	20 14	25 20	59 50	98 95
	CR-UCD	Serif Sans Serif	31 25	50 30	52 44	91 70	100 98

Table 3 :

 3 Employed x-height thresholds

	Character sizes	x-height thresholds (pixels)
	Small	0-30
	Medium	30-55
	Large	55-100