Takwa Chihaoui
email: takwa.chihaoui@esiee.fr

Hejer Jlassi

Rostom Kachouri
email: rostom.kachouri@esiee.fr

Kamel Hamrouni
email: kamel.hamrouni@enit.rnu.tn

Mohamed Akil
email: mohamed.akil@esiee.fr

PERSONAL VERIFICATION SYSTEM BASED ON RETINA AND SURF DESCRIPTORS

Keywords: Biometrics, Retinal verification systems, Optical Disc interest Ring (ODR) method, Speeded up Robust Features (SURF)

Today, Human recognition, especially based on retina, has been an important and attractive topic of scientific research. Most efforts in Biometrics tend to develop more efficient systems which compromise speed and robustness of authentication. In fact, retinal images often suffer from imperfections such as background intensity variation, affine transformations (translation, rotation, scale changes) variations from pattern to other. These defects can seriously affect features extraction in terms of quality and execution time. In this context, in order to overcome these defects, we propose in this paper a novel retinal verification system based on the Speeded Up Robust Features (SURF) extraction. This feature extraction method is so fast and invariant to the affine transformations such as rotation, scale changes and translation. We employ the Optical Disc interest Ring (ODR) method as a preprocessing step in order to further speed up the system and improve the performance. A subset of the VARIA database is used to evaluate the proposed SURF based system. It compromises a high quality with 100% of verification accuracy rate and a time processing very lower than existing verification systems.

I. INTRODUCTION

Automatic Human recognition has been an important and attractive topic of scientific research. Indeed, biometric authentication gets the leading role, since it cannot be borrowed, stolen or forgotten and forging is practically impossible. It is used as the most secure recognition tool. Biometrics is the automatic recognition of a person based on his distinguishing features. In Fact, these are all physical or personal traits automatically measurable, robust and distinctive features that can be used to identify an individual or verify the claimed identity of an individual [START_REF] Polash | Human iris recognition for biometric identification[END_REF]. Various biometric systems have been developed based on fingerprint, palmprint, hand geometry, vein, face, iris, retina and ear. Retinal recognition has received increasing attention in recent years as it provides promising solution to security issues [START_REF] Hill | Retinal identification, in Biometrics: Personal Identification in Networked Society[END_REF] due to its uniqueness, universality, time-invariance and difficulty to be forged. Retinal recognition is a relatively new approach, compared to other biometric features. It includes two modes: Identification which compare each pattern to all others (1:N) and Verification which compare the pattern to others from the same individual (1:1). In this paper, we are interested in the verification mode. Since retinal patterns have highly distinctive traits, the features extracted from retina identify effectively persons, even among genetically identical twins. This pattern will not change through the life of the individual, unless a serious pathology appears in the eye. Most common diseases such as diabetes do not change the pattern in a way that its topology is affected. Hence, this is one of the most secure and valid biometric trait for personal authentication. However, Retinal verification systems can suffer from some problems. Due to the acquisition process of the retinography, retinal images often suffer from imperfections such as background intensity variation, affine transformations (translation, rotation, scale changes) variations from pattern to other. These defects can seriously affect features extraction. In order to overcome these imperfections, we propose, in this paper, a new approach for verifying retina images based on the ODR preprocessing and Speeded Up Robust Features (SURF) description. The verification task in the proposed system is invariant from affine transformations such as rotation, scale changes and translation, and illumination variations. The rest of this paper is organized as follows: Section 2 introduces some existing retinal verification systems. Section 3 presents the proposed system. The pre-processing stage, feature extraction procedure and decision making process is discussed in this Section. Experimental results are presented in Section 4. conclusion is given in the last section.

II. RETINAL VERIFICATION SYSTEMS: RELATED WORKS

Many retina authentication systems have been exploited in order to find the most secure biometric one. Xu et al. [START_REF] Xu | The blood vessel recognition of ocular fundus[END_REF] proposed a new method for recognition. First, they extracted blood vessel skeleton feature using contrast-limited adaptive histogram equalization. After the filtering treatment and extracting shape feature, shape curve of blood vessels was obtained. Shape curve matching was later carried out by means of reference point matching. In their method for recognition, feature matching consisted of finding affine transformation parameters which relates the query image and its best corresponding enrolled image. The computational time of this method is high because a number of parameters should be computed for all possible correspondences between the input and reference images in the dataset. Ortega et al. [START_REF] Ortega | Personal verification based on extraction and characterization of retinal feature points[END_REF] proposed a novel personal verification system based on creases-based method. The feature extraction step consists of generating a Creases image representing the main vessels in order to extract landmarks. After that, they employed a registration process based on the similarity transformation (ST), in order match the reference pattern from the database and the input one. When the biometric patterns were aligned, they measured the degree of similarity by means of a similarity metric to correctly classify authentications in both classes: attacks and clients. Marino et al. [START_REF] Marino | Personal authentication using digital retinal images[END_REF] introduced a novel verification system based on retinal vessel tree pattern. In that work, the whole vascular tree structure was used as the feature pattern. The results showed a high confidence band in the authentication process but with a small database which included only 6 individuals with 2 images for each of them. One of the weak points of this proposed system was the very high computation time of storing and handling a whole image as the biometric pattern. Oinonen et al. [START_REF] Hannu Oinonen | Identity verification based on vessel matching from fundus images[END_REF] proposed a novel method for verification based on minutiae features. The proposed method consisted of three steps: blood vessel segmentation, feature extraction, and feature matching. In practice, vessel segmentation can be viewed as a preprocessing phase for feature extraction. Then, vessel crossings and their orientation information were obtained. These data were matched with the corresponding ones from the comparison image. The method used the vessel direction information for improved matching robustness. The computational time of this method for segmentation, feature extraction, and matching is high. Amin Dehghani et al. [START_REF] Meng | Retinal Identification Based on an Improved Circular Gabor Filter and Scale Invariant Feature Transform[END_REF] proposed a new human recognition (identification and verification) method based on Harris corners and similarity function of retinal images without using vessel segmentation methods. As a preliminary step, they compensated the rotation angle of head or eye movement in front of the retina fundus camera using phase correlation technique. After that, they detected corners in the retina image using Harris corner detector and pruning algorithm. The weakness of this system is that it may be not effective when retinal image are affected by some pathologies. Arathi Arakala et al. [START_REF] Dehghani | Human recognition based on retinal images and using new similarity function[END_REF] proposed an automatic retina verification framework. First, they extracted the vascular pattern of retinal image using a family of matched filters in the frequency domain and morphological operators. Then, they created a spatial graph template based on the extracted pattern. Finally, they compared between two graphs derived from two retinal images using a biometric graph matching (BGM) algorithm which is robust to translation, nonlinear distortion, and small rotations. As a decision step, this system returns a set of distance measures between them in order to distinguish between genuine and imposter comparisons. A common problem on previous approaches is that are variant to affine transformations such as translation, rotation and scale changes, which necessitate a registration phase. This additional phase consists of align the target images with the reference image. Moreover, most of existing approaches suffer from an important execution time. So, in this work, the use of Speeded Up Robust Features (SURF) for feature extraction aims prin-cipally to eliminate the registration phase and consequently to reduce the processing time of the verification system.

III. PROPOSED SURF-BASED RETINAL VERIFICATION

SYSTEM

Depending on several factors, such as the eye location in the objective, patterns may suffer from some deformations. A reliable and efficient model is necessary to deal with these deformations allowing transforming the candidate pattern in order to get a pattern similar to the reference one. The movement of the eye in the image acquisition process basically consists in translation in axis, rotation and sometimes a very small change in scale. In this work, we propose a novel retinal verification system based on the Speeded Up Robust Features (SURF) feature extraction. In the proposed system, segmentation methods that increase computational time were not used. In order to improve the quality and reduce the processing time, a pre-processing step based on the ODR method that enhances input retinal images and extracts an interest ring around the optical disc is needed (Section A). Next, the fast extraction of SURF feature points is explained in Section B. The final step in the verification process is the matching between the keypoints from the reference image and the keypoints from the acquired image based on the test g2NN (Section C).

A. Ring extraction around the optical disc

On the one hand, the presence of noise, the low contrast between vasculature and background, brightness, and the variability of vessel diameter, shape and orientation are the main obstacles in retinal images. On the other hand, SURF characterization is based on the intensity variation between regions, which requires the preservation of the brightest area (the optical disc). In order to overcome these problems, we use our recently proposed method Optical Disc interest Ring (ODR) [START_REF] Mehdi Lajevardi | Retina Verification System Based on Biometric Graph Matching[END_REF] for preprocessing the retinal image in order to improve its quality and extract an interest ring centered at the center of the optical disc as an interest region of retinal image.

B. Feature extraction

In this step, we apply the Speeded Up Robust Features transform (SURF) [START_REF] Takwa | Human Identification System based on the detection of Optical Disc Ring in retinal images[END_REF] to provide a feature description of the retinal image. The output of this step is a 64 elements feature vector for each detected point. SURF is an algorithm for detection of interest areas. It is similar to SIFT descriptor in terms of repeatability and robustness, but exceeds in terms of speed extraction and features vector comparison. Thus, this method is very fast in processing time and is invariant to scaling, rotation and translation. The SURF algorithm involves these steps: an integral representation of the image, a filtering of the image using the Hessian matrix, a suppression of nonmaxima points and a calculation of the descriptors using the Haar wavelet.

1) Integral representation of the image: SURF uses an alternative re-presentation of image called "Integral Image" to speed up the calculation of features. The result is an integral image that has identical dimensions to the input image.

2) Filtering using the Hessian matrix: This step aims on the detection of interest points using the determinant of the Hessian matrix. The Hessian matrix of 2 dimension shown on the equation ( 1) is a square matrix of Gaussian second order derivative with x and y direction of the image f, which is variations of the image around a point.

H(f (x, y)) = df 2 dx 2 df 2 dxy df 2 dyx df 2 dy 2 (1)
The approximated determinant of this Hessian matrix shown below use three filters to calculate the elements of the matrix H(f(x,y)). The use of the Hessian determinant becomes effective with an analysis scales, keeping only the interesting points independently of the observation scale.

3) Suppression of non-maxima: After calculating the Hessian responses SURF search the maxima points which are the interest points of the image on different scales. Then, the nonmaxima points are deleted.

4) Calculation of descriptors :

The last step is to calculate the descriptor vector for each zone of interest. Indeed, for each sub-region 4 * 4, we construct a vector v of 4 elements, defined as in equation ( 2

): v =     dx |dx| dy |dy|     (2)
This vector is a descriptor of the sub-region distribution of gradient. The 16 sub-regions vectors are then aggregated into a single vector of 16 * 4 = 64 which is the descriptor of the region of interest. So, the size of the descriptor vector is 64 ,smaller than SIFT descriptor [START_REF] Mehdi Lajevardi | Retina Verification System Based on Biometric Graph Matching[END_REF] in order to improve the speed comparison . The following figure 3 shows the SURF keypoints distribution on the preprocessed retinal image.

C. Feature matching strategy

This step of matching is performed in the SURF space among the feature vectors of each keypoint to identify similar regions in the sampled image. According to Lowe [START_REF] Bay | Speeded-Up Robust Features (SURF)[END_REF], he limit the classic approach to match keypoints which fixe a threshold on the Euclidean distance between descriptors because of the high-dimensionality of the feature space. However, he used a matching technique called g2NN test to find the best candidate which is based not only on the distance with the first most similar keypoint, but also with the second one; in particular, he used the ratio between the distance to the candidate match and the distance to the second similar feature point. The Principe of this test of matching is that the ratio between the distance of the candidate match di and the distance of the 2nd nearest neighbor di+1 is low in the case of a match (e.g. lower than 0.6) and very high in case of two indiscriminate features (e.g. greater than 0.6). If n is the value in which the procedure stops, each keypoint in correspondence to a distance in d1,. . .,dk where (1 <n< m) is considered as a match for the examined key-point. Finally, by iterating over all keypoints, we can obtain the set of matched points. All the matched keypoints are retained to be linked in the next step. The following figure 4 illustrates an example of SURF keypoints matching between two images.

IV. EXPERIMENTAL RESULTS

All the experiments are implemented in MATLAB, and performed on a PC with a 3.2 GHz CPU and 4 G Byte memory. To evaluate our retinal verification system, the VARIA [START_REF] David | Distinctive Image Features from Scale-Invariant Keypoints[END_REF] which is the largest currently available publicly database, is used for the experiment. It has 233 images from 139 people, with between one to seven samples per person. The images with only one sample and some samples that were identical to each other are excluded from the database. This left 135 retina images from 57 individuals including 42 with 2 samples, 10 with 3 samples, 4 with 4 samples and 1 with 5 samples. The lack of large publicly available retinal databases inhibits large scale and consistent testing of any retinal verification system. The Receiver Operating Characteristic (ROC) curve shows the measured accuracy performance of our biometric system. It illustrates the relationship between (1-FRR) and FAR is shown in Figure5. Our system provides also an AUC (Area Under Curve) value of 100%. Figure 6 shows the False Acceptance Rate (FAR) and the False Reject Rate (FRR) curves. A usual error measure in biometric verification system is the equal error rate (EER) that indicates the error rate where FAR curve and FRR curve are equal. EER is found from these curves as shown in Figure6. In order to evaluate the efficiency of our proposed system in terms of quality, we compare it with some retinal verification systems [3] [4] as shown in Table I.

Fig. 1 .

 1 Fig. 1. Flowchart of our proposed retinal system.

Fig. 2 .

 2 Fig. 2. The result of ODR preprocessing step: (a) Input retinal image and (b) the ring extracted around the optical disc.

Fig. 3 .

 3 Fig. 3. SURF Keypoints distribution: (a) the input interest ring extracted and (b) the ring SURF keypoints characterized.

Fig. 4 .

 4 Fig. 4. Example of g2NN test based SURF keypoints between the clients (authorised) (a) and attacks (unauthorised) (b) authentication cases

Fig. 5 .

 5 Fig. 5. The Receiver Operating Characteristic (ROC) curve of our proposed system.

Fig. 6 .

 6 Fig. 6. The FRR and FAR curve of our SURF-based system.

As shown in Table1, our method verifies 100% of retinal images accurately, as Hannu et al. [START_REF] Hannu Oinonen | Identity verification based on vessel matching from fundus images[END_REF] and Xianjing Meng et al. [START_REF] Meng | Retinal Identification Based on an Improved Circular Gabor Filter and Scale Invariant Feature Transform[END_REF] systems . This is due to the invariance to the rotation and translation ensured by SURF. As shown in table I, the Equal Error Rate (EER) of our system achieved the 0% however other systems [3] [4] suffer respectively from an error rate (Equal Error Rate) 0.19% and 0.01%. The most gain of our proposed system compared to other verification systems is in term of processing time. In fact, as shown in table II, our system is the fastest one with an average of processing time of 0.03s, compared to respectively Xu et al. [START_REF] Xu | The blood vessel recognition of ocular fundus[END_REF] system which is very consuming time with approximately 45s and Xianjing Meng et al. [START_REF] Meng | Retinal Identification Based on an Improved Circular Gabor Filter and Scale Invariant Feature Transform[END_REF] with 7.87s. The huge execution time of the system [START_REF] Xu | The blood vessel recognition of ocular fundus[END_REF] is due to the complicate matching process which compute a number of rigid motion parameters should be computed for all possible correspondences between every the input and enrolled images.

V. CONCLUSION

In this paper, we present an automatic retina verification system based on Surfing process. It is described as follows. Firstly, the pre-processing phase consists of the extraction of a disc interest region around the optical disc to be an input for the second phase. Secondly, the features of the extracted region are extracted by the SURF descriptor. Finally, the test g2NN is applied to compute the number of matched keypoint pairs and classify identical retinal images. Our experimental results show that this work compromises an exact personal verification and a very lower computational complexity. The accuracy of the proposed system for human verification is 100%.


Review & Perspective for Distance Based Clustering of Vehicle Trajectories

Philippe C. Besse, Brendan Guillouet, Jean-Michel Loubes, and Franc ¸ois Royer, Abstract-In this paper we tackle the issue of clustering trajectories of geolocalized observations based on distance between trajectories. We first provide a comprehensive review of the different distances used in the literature to compare trajectories. Then based on the limitations of these methods, we introduce a new distance: Symmetrized Segment-Path Distance (SSPD). We compare this new distance to the others according to their corresponding clustering results obtained using both the hierarchical clustering and affinity propagation methods. We finally present a python package : trajectory distance, which contains the methods for calculating the SSPD distance and the other distances reviewed in this paper.

Index Terms-Trajectory clustering

I. INTRODUCTION

A TRAJECTORY is a set of positional information for a moving object, ordered by time. This kind of multidimensional data is prevalent in many fields and applications, for example, for understanding migration patterns through studying trajectories of animals, predicting meteorology with hurricane data, improving athletes performance, etc. Our study concentrates on vehicle trajectories within a road network. The growing use of GPS receivers and WIFI-embedded mobile devices equipped with hardware for storing data enables the collection of an enormous amount of data that can be used to extract relevant information in order, for instance, to find the optimal path to go from point A to point B, detect abnormal behavior, optimize the traffic flow in a city, predict the next location or final destination of a moving object etc. This aims actually to build from the data the different features that characterize the different daily movement of the vehicles on the road network. For this, we consider clustering methods for trajectories. Clustering techniques aim to regroup similar trajectories together into groups that are different from one another. The complexity of trajectory makes this a challenging task as objects can move along many different paths in a given area, moreover the road network of a highway does not have the same complexity as that of a city, and finally the road network in a city differs between the suburbs and downtown. In addition, the speed of an object varies between regions, and between paths taken within a single region. Even within the same path the speed depends on exogenous variables such as the time of day, or whether it is a weekday or the weekend.

Several methods can be used to cluster trajectories. We will focus on distance based trajectory clustering but other specific methodologies have been investigated. Dealing with the functional properties of trajectories considered as continuous function of time Gaffney (2009, [START_REF] Gaffney | Trajectory clustering with mixtures of regression models[END_REF]), Vasquez et al. (2004[2]), Hu et al. (2006[3]) successfully apply trajectory clustering methods on video-stream trajectories. Gariel et al. [START_REF] Gariel | Trajectory clustering and an application to airspace monitoring[END_REF] also use the continuous definition of trajectory to re-sample the trajectories and obtain time series of equal length. A principal components analysis is then applied on these new trajectories to obtain principal components and finally cluster them. This method is applied to airplane routes. All these methods take into account both the spatial and the temporal aspect of the trajectory, they are not adapted to the vehicle trajectories constrained to a road network whose time progression is very irregular. Rinzivillo et al. (2008,[5]) and Kim et al. [START_REF] Kim | Spatial and temporal characterization of travel patterns in a traffic network using vehicle trajectories[END_REF]) propose density-based methods. Both require the definition of a density parameter and a minimum cluster size which implies an extensive knowledge of the studied area or a precise question to obtain good results. Hence, these methods are hard to automatically adapt from one dataset to another. Finally Lee et al. (2007,[7]) and Wu et al. (2014[8]) propose to use clustering methods on trajectory line segments to enable the detection of important areas of flow, though this does not consider the trajectory as a whole path. Our main objective here is to detect the main path and traffic flow which can later be used to study the different behaviors along these paths.

In this context, the goal of this work is to construct, in a data driven way, a collection of trajectories that model the behaviors of car drivers. These models are learned from a data set of car locations. In this work we focus on clustering trajectories having similar paths. This clustering is based on the comparison between trajectory objects, and as such a new definition of distance between the studied trajectory objects is required.

A large amount of work has been done to give new definitions of trajectory distance. Tiakas et al.(2009[9]) , Rossi et al. (2012 [START_REF] El Mahrsi | Graph-based approaches to clustering network-constrained trajectory data[END_REF]), Han et al.(2015[11]) or Hwang et al.(2005[12]) propose road network based distances. They assume that the trajectories studied are perfectly mapped on the road network. However, this task is strongly dependent on the precision of the GPS device. When the time interval between two GPS locations is significant, several paths on the graph are possible between locations, especially when the network is dense. Moreover it requires the knowledge of the road network. Here, we focus on entirely data driven methods without any a priori information. Several methods have been used to cluster data set of trajectories. Clustering methods using Euclidean distance lead to inaccurate results mainly because trajectories have different lengths. Hence, several methods based on warping distance have been defined , Berndt (1994[13]), Vlachos et al. (2002[14]), Chen et al. (2004 [START_REF] Chen | On the marriage of lp-norms and edit distance[END_REF]), and Chen et al. (2005 [START_REF] Chen | Robust and fast similarity search for moving object trajectories[END_REF]). These methods reorganize the time index of trajectories to obtain a perfect match between them. Another approach is to focus on the geometry of the trajectories, in particular their shape. Shape distances like Hausdorff and Fréchet distances can be adapted to trajectories but fail to compare them as a single entity. Lin et al. (2005[17]) proposed a method based exclusively on the shape of the trajectory but at high computational cost.

In section II the papers definitions, notations and problem statement are introduced. In section III several distances on trajectory are studied and compared. A new distance will be presented in section IV: the Symetrized Segment-Path Distance (SSPD). SSPD is a shape-based distance that does not take into account the time index of the trajectory. It compares trajectories as a whole, and is less affected by incidental variation between trajectories. It also takes into account the total length, the variation and the physical distance between two trajectories. For all these different distances, we obtain different clusterings. So we can compare the distance on these results. The choice of clustering used is detailed section V. The SSPD and the other studied distances were implemented in a python package, trajectory distance available on github. The presentation of this package and the experimental evaluation of these distances with the chosen clustering techniques on some trajectory sets are analyzed in section VI.

II. MODEL FOR TRAJECTORY CLUSTERING A. Trajectory

A continuous trajectory is a function which gives the location of a moving object as a continuous function of time.

In our case we will only consider discrete trajectories defined here after. Definition 1. A trajectory T is defined as

T : ((p 1 , t 1 ), . . . , (p n , t n )), where p k ∈ R 2 , t k ∈ R ∀k ∈ [1 . . . n],
∀n ∈ N and n is the length of the trajectory T .

The exact locations between time t i and t i+1 are unknown. When these locations are required, a piece wise linear representation is used between each successive location p i and p i+1 resulting in a line segment s i between these two points. This new representation is called a piece wise linear trajectory. In this representation, no assumption is made about time indexing of segment s i . Definition 2. A piece wise linear trajectory is defined as T pl : ((s 1 ), . . . , (s n-1 )) , where s k ∈ R 4 and n pl is the length of the trajectory.

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

n pl = i∈[1...n-1] p i p i+1 2 .
The notation used in this paper are summarized in Table I.

B. Distance

There are many ways to define how close two objects are from one another. Beyond the notion of mathematical distance, many functions can be used to qualify this dissimilarity. The terminology used in literature to define them is not completely standardized. Therefore we will use the definition established in Deza et al. (2009[18]) as a reference. Definition 3. Let T be a set of trajectories. A function d : T ×T → R is called a dissimilarity on T if for all T 1 , T 2 ∈ T :

• d(T 1 , T 2 ) ≥ 0 • d(T 1 , T 2 ) = d(T 2 , T 1 ) • d(T 1 , T 1 ) = 0
If all of these conditions are satisfied and d(T 1 , T 2 ) = 0 =⇒ T 1 = T 2 d is considered to be a symmetric. If the triangle inequality is also satisfied, d is called a metric. These notations are summarized in Table II.

X indicates the required properties for each distances, while * indicates properties that are automatically satisfied (by the presence of the other required properties for the metric).

C. Desired properties of clustering and distances

Our aim is to regroup trajectories sharing similar behaviour. We want that trajectories in the same cluster, take similar paths. Hence our goal is to define a clustering method that will regroup trajectories • with similar shape and length • which are physically close to each other • which are similar as a whole with more than just similar sub-parts • all of these properties should be considered without regard to their time indexing Moreover we want to design a very general procedure which is able to treat all trajectories data, without prior knowledge of the particular geographical location where they are collected. To obtain such clustering, the goal of this work is to find a distance that respects such properties and to succeed in extracting these features. Actually, the desired distance should have the following properties,

• it compares trajectories as a whole • the compared trajectories can be of different lengths,

• the time indexing can be very different from one trajectory to another • the trajectories can have similar shapes but can be physically far from each other and vice versa

D(T 1 , T 2 ) ≥ 0 X X * Symmetry D(T 1 , T 2 ) = D(T 2 , T 1 ) X X X Reflexivity D(T 1 , T 1 ) = 0 X * * Triangle Inequality D(T 1 , T 3 ) ≤ D(T 1 , T 2 ) + D(T 2 , T 3 ) X Identity of indiscernible D(T 1 , T 2 ) = 0 =⇒ T 1 = T 2 X X
• extra parameters should not be required.

III. DISTANCE ON TRAJECTORIES: A REVIEW

Three main kind of distances have been introduced in the literature. The first uses the underlying road network, Network-Constrained Distance. These distances will not be detailed in this paper. They assume that the road network is known and that trajectory data are perfectly mapped on it. Distances that do not use the underlying road network can also be classified into two categories: those who only compare the shape of the trajectory, Shape-Based Distance and those who take into account the temporal dimension, Warping based Distance.

Performance of clustering algorithms using these distances will be compared in section II, as well as their computation cost and their metric properties.

A. Warping based Distance

Euclidean distance, Manhattan distance or other L p -norm distances are the most obvious and the most often used distances. They compare discrete objects of the same length. They can be used to look for common sub-trajectories of a given length but they cannot be used to compare entire trajectories. Moreover, these distances will compare locations with common indexes one by one. At a given index i, location p 1 i of trajectory T 1 will be compared only to location p 2 i of trajectory T 2 . However, these locations can be strongly different according to the speeds of the trajectories. Hence, it makes no sense to compare them without taking this into account. This problem is also common in time series analysis and not restricted to trajectory analysis. Warping distance aims to solve this problem. To accomplish this, they enable matching locations from different trajectories with different indexes. Then, they find an optimal alignment between two trajectories, according to a given cost δ between matched location. Several warping based distances have been defined. DT W (Berndt et al., (1994 [13])) and later LCSS (Vlachos et al., 2002[14]), EDR (Chen et al., 2005[16]) and ERP (Chen et al., 2005[15]). These distances are defined the same way, but they use different cost functions.

In order to define a warping distance, two compared time series trajectories, T i , T j , are arranged to form a n i × n j grid G. The grid cell, g k,l , corresponds to the pair (p i k ,p j l ). Definition 4. A warping path, W = w 1 , . . . , w |W | , crosses the grid G such that

• w 1 = g 1,1 , • w |W | = g n i ,n j , • if w k = g ki,kj , then w k+1 is equal to g ki+1,kj , g ki,kj +1
or g ki+1,kj +1 .

The order of the locations in a trajectory are maintained but they can be repeated, deleted or replaced by an arbitrary value, a gap, along the warping path. The distance is then computed by minimizing or maximizing the sum of a given cost δ between all pair of locations that make a warping path W , for all existing warping paths. Definition 5. A warping distance is defined as

D(T i , T j ) = min W |W | k=1 δ(w k ) , or = max W |W | k=1 δ(w k ) , (1)
where δ(w k ) = δ(g ki,kj ) = δ(p i ki , p j kj ), is the cost function and W is a warping path.

They are generally computed by dynamic programming. Table III displays the cost functions as well as the dynamic formulation of these distances.

Contrary to the three other distances, LCSS is a similarity. The exact similarity used in Vlachos et al., 2002[14] is S(T i , T j ) = LCSS(T i ,T j ) min{n i ,n j } , which is between 0 and 1. We will then use the distance

DLCSS(T i , T j ) = 1 -S(T i , T j ),
to compare distances to each other.

The metric types of these distance functions, and computational cost for the four methods are summarized in table IV.

1) Comparisons:

• All of these distances handle local time shifting.

• The cost function δ uses the Euclidean distance. Some of these distances have been defined using a L1-norm, but Euclidean distance is more adapted for real values. • LCSS and EDR's cost function count the number of occurrences where the Euclidean distance between matched location does not match a spatial threshold, ε d . The former counts similar locations, the latter the difference. This threshold makes the distance robust to noise. However, it has a strong influence on the final results. If the threshold is large, all the distances will be considered similar and if low, only those having very close locations will be considered similar. • In comparison, ERP and DT W add a weighting to these differences by computing the real distance between

N AM E (p1, p2) = N AM E(T i , T j ) = DTW p1p2 2 =                0 if n i = n j = 0 ∞ if n i = 0 or n j = 0 δ DT W (p i 1 , p j 1 )+ min DT W (rest(T i ), rest(T j )), DT W (rest(T i ), T j )), DT W (T i , rest(T j ) otherwise NAME LCSS    1 if p1p2 2 < ε d 0 if p1 or p2 is a gap 0 otherwise =          0 if n i = 0 or n j = 0 LCSS(rest(T i ), rest(T j )) + δ LCSS (p i 1 , p j 1 ) if δ LCSS (p i 1 , p j 1 ) = 1 max LCSS(rest(T i ), T j )) + δ LCSS (p i 1 , gap), LCSS(T i , rest(T j )) + δ LCSS (gap, p j 1 ) otherwise EDR    0 if p1p2 2 < ε d 1 if p1 or p2 is a gap 1 otherwise =                n i if n j = 0 n j if n i = 0 EDR(rest(T i ), rest(T j )) if δ EDR (p i 1 , p j 1 ) = 0 min EDR(rest(T i ), rest(T j )) + δ EDR (p i 1 , p j 1 ), EDR(rest(T i ), T j )) + δ EDR (p i 1 , gap), EDR(T i , rest(T j ) + δ EDR (gap, p j 1 ) otherwise ERP    p1p2 2 if p1, p2 are not gaps p1g 2 if p2 is a gap gp2 2 if p1 is a gap =              n i k=1 p i k g 2 if n j = 0 n j l=1 p j l g 2 if n i = 0 min ERP (rest(T i ), rest(T j )) + δ ERP (p i 1 , p j 1 ), ERP (rest(T i ), T j )) + δ ERP (p i 1 , gap), ERP (T i , rest(T j ) + δ ERP (gap, p j 1 )
otherwise

LCSS distance O(n 2 ) EDR symmetric O(n 2 ) ERP metric O(n 2 )
the locations. In this sense they can be viewed as more accurate.

• ERP is the only distance which is a metric regardless of the L p norm used, yet it works better for normalized sequences, especially for defining the gap value g. It does not apply for vehicle trajectories. • In addition, these distances may include a time threshold, ε t . Thus, two locations will not be compared if the difference between their time indexing is too large. However, it is very hard to estimate the value of this threshold when comparing trajectories due to the presence of noise. 2) Pros and Cons: The main advantage of these distances is that they enable comparison of sequences of different lengths.

The two main limitations of warping based distance are the following • Warping methods are based on one-to-one comparison between sequences. Hence, it often requires the choice of a particular series that will be used as a reference, onto which all other sequences will be mapped. The index of two sequences being compared should be well balanced in order to best capture the variability, for instance, in order to detect if there were accelerations and decelerations during the measurement of the time series. Hence the choice of the reference sequence is very important. • The performance of the usual methods based on warping techniques is hampered by the large amount of noise inherent to road traffic data, which is not the case when examining time series. Instead of correcting the time index, the solution is to use distances that have the effect of time removed.

B. Shape-Based Distance

These distances try to catch geometric features of the trajectories, in particular, their shape. Among Shape-Based Distances, the Hausdorff distance (Hausdorff, 1914 [19]), and the Fréchet distance (Fréchet, 1906[20]) are likely the most well known.

1) Hausdorff:

The Hausdorff distance is a metric. It measures the distance between two sets of metric spaces. Informally, for every point of set 1, the infimum distance from this point to any other point in set 2 is computed. The supremum of all these distances defines the Hausdorff distance. Definition 6. The Hausdorff distance between two sets of metric spaces is defined as

Haus(X, Y ) = max{sup x∈X inf y∈Y xy 2 , sup y∈Y inf x∈X xy 2 }.
This distance is complicated and resource intensive to compute when applied to most existing sets. But in the case of polygonal curves like trajectories, some simplification can be made due to the monotonic properties of a segment. Distance from a point p to a segment s is defined as follows.

Definition 7. 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 .

Hence, the Hausdorff distance between two line segments is

D Haussdorf (s 1 i1 , s 2 i2 ) = max{ sup p∈s 1 i 1 D ps (p, s 2 i2 ), sup p∈s 2 i 2 D ps (p, s 1 i1 )} = max{ D ps (p 1 i1 , s 2 i2 ), D ps (p 1 i1+1 , s 2 i2 ), D ps (p 2 i2 , s 1 i1 ), D ps (p 2 i2+1 , s 1 i1 )}.
Indeed, a segment is monotonic. As seen in Fig. 1, the supremum of the P oint -to -Segments distance from any point of a segment s 1 i1 to a segment s 2 i2 occurs at one of the end points of the segment s 1 i1 . The Hausdorff distance between two trajectories can then be computed with the following formula. Definition 8. Hausdorff distance between two discrete trajectories.

D Haussdorf (T 1 , T 2 ) = max max i 1 ∈[1...n 1 ] j 2 ∈[1...n 2 -1] {Dps(p 1 i1 , s 2 j2 }, max j 1 ∈[1...n 1 -1] i 2 ∈[1...n 2 ]
{Dps(p 2 i2 , s 1 j1 } .

The Hausdorff distance can then be computed in a O(n 2 ) computational time.

2) Frechet and discrete Fréchet : The Frechet distance measures similarity between curves. It is often known as the "walking-dog distance". Imagine a dog and its owner walking on two separate paths without backtracking from one endpoint to one other. The Fréchet distance is the minimum length of leash required to connect a dog and its owner. While the Hausdorff distance takes distance between arbitrary points, the Fréchet metric takes the form of the two curves into account.

Definition 9. The Fréchet distance between two curves is defined as

D F rechet (A, B) = inf α,β∈X max t∈[0,1] A α(t) , B β(t) 2 .
Similar to the Hausdorff distance, the Frechet distance is a metric. It is also resource intensive. Alt et al. (1995[21]) developed an algorithm measuring the exact Fréchet distance for polygonal curves based on the free space definition.

Definition 10. The free space F (T 1 , T 2 ) between two trajectories is the set of all pairs of points whose distance is at most .

F (T 1 , T 2 ) := {(p 1 , p 2 ) ∈ (T 1 , T 2 )}| p 1 , p 2 2 ≤ }.
The Fréchet distance between two trajectories T 1 and T 2 is the minimum value of for which a curve exists within the corresponding F from (p 1 0 , p 2 0 ) to (p 1 n 1 , p 2 n 2 ) with the property of being monotone existing in both trajectories. Computing the Fréchet distance means finding the minimum value of . By exploiting the monotonic property of the segments and the definition of free space, this task can be accomplished more efficiently.

Indeed, the Frechet distance between segments is equal to the Hausdorff distance between segments, i.e.

D F rechet (s 1 i1 , s 2 i2 ) = max{ D ps (p 1 i1 , s 2 i2 ), D ps (p 1 i1+1 , s 2 i2 ), D ps (p 2 i2 , s 1 i1 ), D ps (p 2 i2+1 , s 1 i1 )} = i1,i2 .
To compute the Frechet distance between trajectories T 1 and T 2 , we need only look among the set E of Frechet distances between all pairs of segments of T 1 and T 2 .

E = { i1,i2 for (i 1 , i 2 ) ∈ ([1 . . . n 1 -1] × [1 . . . n 2 -1])}.
This simplification enables us to compute the Frechet distance between trajectories T 1 and T 2 in O(n 2 log(n 2 )). We highlight that this computational cost is higher than all the other calculation methods of the studied distances.

Eiter et al. (1994 [START_REF] Eiter | Computing discrete fréchet distance[END_REF]) describes an approximation of this distance for polygonal curves called the discrete Fréchet distance. This distance is close to the definition of the warping based distance.

Definition 11. The discrete Fréchet distance is defined as

D F rechet Discr ((T 1 , T 2 ) = min W { max k∈[1...|W |] w k 2 }.
with W being the warping path defined in definition 5. The discrete Fréchet distance can be computed in O(n 2 ) time.

This distance is bounded as follows.

Theorem 1. For any trajectories T i and T j [22]

D F rechet (T i , T j ) ≤ D F rechet Discr ((T i , T j ) ≤ D F rechet (T i , T j )+ Where, = max{ max k∈[1...n i -1] { p i k p i k+1 2 }, max l∈[1...n j -1]
{ p j l p j l+1 2 }}.

3) One Way Distance: Lin et al. 2005 [START_REF] Lin | Shapes based trajectory queries for moving objects[END_REF] defines the One-Way-Distance, OWD, from a trajectory T i to another trajectory T j as the integral of the distance from points of T i pl to trajectory T j pl divided by the length of T i pl Definition 12. The OWD distance is defined as

D OW D (T i , T j ) = 1 n i pl p i ∈T i pl D point (p i T j )dp i ,
where D point (p, T ) is the distance from the point p to the trajectory T so that

D point (p, T ) = min q∈T pl pq 2 . The OWD distance is not symmetric, but D SOW D (T i , T j ) = (D OW D (T i , T j ) + D OW D (T j , T i ))/2
is. This distance is a symmetric because it does not satisfy the triangle inequality.

Lin et al. [START_REF] Lin | Shapes based trajectory queries for moving objects[END_REF], have defined two algorithms to compute the OWD in case of piecewise linear trajectories.

• The first consists of finding the parametrized OWD function D OW D (s i k , T j ) from a segment s i k of T i pl to all segments s j of T j pl and for all segments of T i pl

D OW D (T i , T j ) = 1 n i pl n i -1 k=1 D OW D (s i k , T j ). p i k p i k+1 ,
with a complexity of O(n 2 log(n)).

• The second uses a grid representation of the trajectory.

As we see in Fig. 2, the space is discrete Trajectories are defined as the succession of grids they have crossed.

Definition 13. A grid representation trajectory is defined as T grid := (g 0 , . . . , g n grid ),

where g n are cells of the discrete space. [START_REF] Lin | Shapes based trajectory queries for moving objects[END_REF]) is by far the one that best meets our requirements. It compares trajectories as a whole, taking into account their shapes and their physical distances, the required features for our distance. However, its complexity makes it computationally slow. The algorithm for grid representation is faster. Its computational time is O(mn). Yet it does not take into account the computation time required for matching the trajectory to the grid. Moreover, the size of the grid chosen strongly influences the final result and makes it imprecise. Furthermore, the distance gives the same "weight" to all points defining the trajectory: points directly issued from the GPS location, and points which compose the piece wise linear representation. The greater the length of the segment s is, the stronger its influence on the trajectory is. The more separated the endpoints of a segment s are, the less confident the interpolation between them is.

In the following section, a new distance will be established inspired from both the OWD and the Hausdorff distances.

IV. A NEW DISTANCE : SYMMETRIZED SEGMENT-PATH DISTANCE (SSPD) In this section, we define a new shape based distance, the Symmetrized Segment-Path Distance, and we compare it to other shape based distances. We propose SSPD in order to fulfill the desired properties defined in section II-C.

Like the Hausdorff distance (Definition 8), the definition of SSPD is based on the P oint -to -Segment distance (Definition 7). From this definition, we define the P oint-to-T rajectory distance, D pt , from a point p to a trajectory T like the minimum of distances between this point and all segments s that compose T (Figure 4). The Segment-Path distance from trajectory T 1 to trajectory T 2 is the mean of all distances from points composing T 1 to the trajectory T 2 (Figure 5). Definition 14. SP D distance is defined as

D SP D (T 1 , T 2 ) = 1 n 1 n1 i1=1 D pt (p 1 i1 , T 2 ).
where,

D pt (p 1 i1 , T 2 ) = min i2∈[0,...,n2-1] D ps (p 1 i1 , s 2 i2 ).
D SSP D (T 1 , T 2 ) = D SP D (T 1 , T 2 ) + D SP D (T 2 , T 1 ) 2 .
In definitions 14 and 15, distances SPD and SSPD are computed by taking the mean of the Point-to-Trajectory distance and the SPD distance. If the maximum is used instead of the mean, one recovers the Hausdorff function between two trajectories. Computing only one distance between two locations makes it very sensitive to noise. Yet our method computes the mean of such quantities which makes it less so. For example, for the trajectories in Fig. 3, the SSPD distance between T 1 and T 2 is smaller than the distance between T 1 and T 3 or T 2 and T 3 (D(T 1 , T 2 ) = 0.58, D(T 1 , T 3 ) = 1.5, D(T 2 , T 3 ) = 2.03).

Proposition 2. SSDP is a symmetric.

Proof. SSDP is a sum of Euclidean distances. By definition SSDP is greater or equal to 0. By definition 15, SSDP is symmetric. Finally theorem 1 states that, if D SDP (T 1 , T 2 ) = 0, T 1 is a sub trajectory of T 2 . Therefore if D SSDP (T 1 , T 2 ) = 0, both D SDP (T 1 , T 2 ) = 0 and D SDP (T 1 , T 2 ) = 0, and

T 1 = T 2 . SSDP is then a symmetric.
SSDP is quite similar to OWD but its definition resolves most of the problems of OWD regarding the desired properties defined in II-C

• The points coming from the interpolation of two observed locations of a trajectory are less trustworthy that the real observations. Hence, it is natural to give more weight to the observed points. • SSPD distance does not require any additional parameters such as a threshold or a grid to be computed. • Its computation cost is O(n 2 ). It only depends on the number of locations.

V. CLUSTERING

To evaluate these different distances, we will study different clustering methods obtained with the same algorithm but with distances computed using all previous distances. The different selected clustering methods and the quality of cluster criterion are examined in this section.

A. Methods

The choice of the clustering method is restricted by the characteristics of the trajectory object. Indeed, trajectories have different lengths which complicates an easy definition of a mean trajectory object. The k-means method cannot be used on our trajectory set, nor spectral clustering method. k-medoid can be used but an efficient algorithm, like partitioning around medoids, or dbscan method, require a valid metrics. Indeed, these algorithms are based on nearest neighbor and require the distance used to be known in order to satisfy the triangular inequality. Most of the studied distances, SSPD, LCSS, DTW, are not metrics. In this way, dbscan or partitioning around medoids algorithms will not be used. Moreover, dbscan depends on two extra parameters that are hard to estimate in this case.

To perform the clustering of the trajectories, we will focus on two methodologies : hierarchical cluster analysis (HCA) and affinity propagation (AP). As a matter of fact, HCA and AP can use distance/similarity which does not satisfy the triangle inequality. We point out that the choice of the clustering method is restricted to the trajectory object we deal with. Actually, trajectories have different lengths. HCA and AP are both methods which only require the distance/similarity matrix, and thus can cluster objects of different lengths. Both of these methods will be used to evaluate our distance.

B. Quality criterion of cluster result

A clustering algorithm aims to gather objects into homogeneous groups that are far one from another. Hence, the quality of a clustering is usually evaluated by looking at the between and within variance of the obtained clusters. On the one hand, the within variance shows the spread of elements belonging to the same groups. Because we want the elements of the same groups to be as similar as possible, we want the within variance to be as small as possible. On the other hand, we want objects that belong to different groups to be as far as possible from one another. Hence, the between variance, which shows the spread between clusters, should be as big as possible. The definition of within group variance and between group variance require the definition of a mean object. In our case, they cannot be computed here because of the impossibility of computing the mean of the trajectory object. Therefore, we approximate this mean by considering an exemplar, T ex , of a set of a trajectory T of length n T , defined as:

T ex T = min T i i∈[0...n T ] n T j=1 j =i D(T i , T j ) .
By using the exemplar definition, we can define new criteria to approximate the between and within variance : the Between-Like and the Within-Like Criteria. Let C 1 , . . . , C K be a set of clusters of T , the Between-Like and the Within-Like criteria are defined as: Definition 16. Between-Like and Within-Like

BC = K k=1 D(T ex T , T ex C k ), W C = K k=1 1 |C k | T i ∈C k D(T ex C k , T i ).
The Within-Like criterion shows the spread of elements belonging to the same cluster while the Between-Like criterion shows the spread between clusters. As for the variance, for a given number of clusters, we want the Within-Like criterion to be as small as possible, and the Between-Like criterion to be as big as possible.

These quality criteria can also be used to select the number of clusters. Indeed, the Within-Like criterion decreases with the number of clusters, while the Between-Like criterion increases with the number of clusters. Hence, we are looking for a trade off between these criteria and the number of clusters. We choose the number of clusters so that adding one more cluster does not decrease significantly the Within-Like criteria.

VI. EXPERIMENTAL EVALUATION

In this section, we evaluate and compare 7 distances LCSS, DTW, Hausdorff, Fréchet, Discrete Fréchet, OWD grid and the SSDP. We present the python package trajectory distance where the distances have been implemented. We compare their computational cost and the results of the application of different clustering techniques for each of them. We also use python for the implementation of the chosen clustering algorithms, the sklearn library for affinity propagation and scipy library for hierarchical clustering analysis. For the latter, weighted, average, ward, complete and single linkage criteria are compared.

A. A python package : trajectory distance

All distances have been implemented in python and are available in the trajectory distance package available on github from this url : https://github.com/bguillouet/traj-dist. Cython is a language which enables the declaration of static variables and to use the C math library. All distances but OWD grid are based on Euclidean distance for point to point computation. To compute OWD grid, we implemented the grid algorithm defined in [START_REF] Lin | Shapes based trajectory queries for moving objects[END_REF]. The trajectories need to be mapped in a grid space. We use the Geohash system to accomplish this task. It is based on a hash function which subdivides the geographical space into a grid. Different precision parameters of this grid are available from 1 (5009.4km × 4992.6km cell) to 12 (3.7cm × 1.9cm). The hash function enables quick mapping of latitude/longitude coordinates to the appropriate grid. Once the trajectory locations have been mapped to a grid, we implemented a algorithm to find all grids they crossed between two locations and to obtain the grid representation of the trajectory as defined in Definition 13.

B. The Data

The data we used are GPS data from 536 San-Francisco taxis over a 24-day period. These data are public and can be found on [START_REF] Piorkowski | CRAW-DAD data set epfl/mobility (v. 2009-02-24)[END_REF]. We extracted a subset of this data as shown Fig. 6.

This subset is a blend of 2574 trajectories. All have the same pickup location, the Caltrain station, and all have a drop-off location in downtown San-Francisco.

C. Computation cost

In Table VI we can observe the computation time needed to compute the matrix distance for a subset of 100 taxis trajectories of the studied subset described in section VI-B. Trajectories are composed of 3 to 39 locations, most having around 10.

Fréchet distance is the distance that requires the most computation time. It is the only method that runs in O(n 2 log(n 2 )). DTW, LCSS and Discrete Fréchet distances are the fastest computed methods, all having the same order of computation time. This three methods require computing the Euclidean distance between each pair of points that compose the two trajectories. They only differ by their cost function. As explained in section IV, Hausdorff and SSPD distance also are computed in the same way. Hausdorff uses the maximum of the Point-to-Trajectory distance and SSPD the mean, which explains why they have almost the same computation time. If they both have the same complexity as that of DTW, LCSS and Discrete Fréchet, log(n 2 ), they require more operations, which explains their relative slowness. The time computation of the OWD grid is strongly dependent on the precision we choose. Indeed, a precision of 5 resolves the studied space into a 3 × 3 grid space, while precisions of 6 and 7 resolve a grid space of 10 × 5 and 34 × 25 respectively. This implies that the number of cells required to represent the trajectory is very different from one precision to another. With precision 5, 1 to 3 cells are needed to represent the trajectory. With precision 7, the shortest trajectory is represented with only 1 cell but the longest needs 47 cells to fully describe it. Therefore, with precision 5, OWD grid is faster than SSPD, but with precision 6 and 7, OWD grid is 3 times and more than 20 times faster than SSPD distance. A good trade-off needed to be found to have a precision which enabled to adequately describe the trajectories in the new grid space, within a reasonable time computation.

D. Analysis of the clustering method

In Fig. 7 we observe the evolution of the Within-and the Between-like criteria described in section V for the distance SSPD and for the selected methods AP and HCA. Both the Between-Like and the Within-Like criteria are shown because the sum of these two criteria is not constant as opposed to the sum of the between and within variance.

The HCA single method gives much worse results than the other methods, regardless of the number of cluster. All other HCA methods display the same evolution of the studied Fig. 7: Evolution of the Between-Like (a) and Within-Like (b) criteria depending on cluster size for clusterings obtained with SSPD. On (c), Within-Like criterion is displayed with a reduced scale of 0 to 50. criteria with respect to the cluster size. A plateau can be observed starting at cluster sizes between 10 and 20. Adding more clusters does not decrease significantly the Within-Like Criterion. We can see Figure 7-c, that the HCA ward give better results than the other HCA method for all the different number of clusters. The same conclusions can be made regardless of the distance used.

AP give the best results. However the latter does not achieve clustering using any less than 36 clusters, a potentially large number, given that HCA ward achieves clustering around 15, and with good Within-Like criterion results. Moreover, the minimum cluster size found by the AP method differs significantly according to the distance used. No less than 21 clusters are found with the DTW distance and 54 with Hausdorff. This is the main inconvenient of this method.

The HCA Ward method and the AP method with the preference parameter fixed to the minimum of the computed matrix distance will be used to compare the studied distances in more details.

E. Analysis of the distances

We can observe the evolution of the Within-Like and the Between-Like criteria for the two selected clustering methods as well as for all studied distances. The HCA WARD results are display in Fig. 8, and the AP results in Fig. 9.

In Figure 8, the evolutions of the two criteria when the number of cluster increase is similar for all the distances but LCSS. For instance, the curves which represent the Within-Like criterion decrease quickly when adding more cluster for a low number of cluster. Then they all reach a point when adding more clusters does not decrease significantly the Within-Like criterion. Theses point vary from one distance to another. It is However, the same conclusion can be made in terms of Within-Like and Between-Like criteria regardless of the number of cluster.

The Warping-based distances, LCSS and DTW, give the poorest results with LCSS being significantly worse than DTW. The two shape-based distances Frechet and Hausdorff give better results. The evolution of their criteria is very similar to one another. The Discrete Fréchet distance is between these two types of distances. These results confirm that shape-based distances are better adapted than warping-based distances for our objectives.

OWD grid gives the best results. It has the lowest value of Within-Like Criterion for all cluster sizes using both HCA WARD and AP clustering methods. The results for precision 5 are not displayed here. The discretization of the space is too inaccurate and did not provide good clustering. Precision 7 gives slightly better results than precision 6. This shows that increasing the precision parameter yields better results. However, when there are more than 15 clusters found, the within and Between-Like criteria are almost the same. In section VI-C we have seen that the computation time to compute the distance with precision 7 is seven times higher than the computation time with precision 6. Hence, we need to look for the optimal criteria to find a good trade off between good clustering results, and reasonable computational time. The choice of this criteria is a strong disadvantage of the OWD method, because it implies to look for the best precision parameter for each data set.

Finally, the new distance SSDP is the distance which best approaches the results found with OWD grid, regardless of the number of cluster. But unlike with OWD grid, we do not need to look for the optimal precision parameter, in order to compute it, nor to map the trajectory to a new space. This enables our distance to be more easily adapted to different subsets of trajectories.

We observe the visual results for this distance and both AP and HCA ward clustering methods, in Fig. 10, and the isolated clusters, in Fig. 11. For the HCA ward method, we display the clustering result obtained with 15 clusters because we have seen Section VI-D that a plateau can be observed on the evolution of the Within-Like criteria with respect to the number of cluster starting at cluster sizes between 10 and 20 for the SSPD. For the AP, clustering results with 36 clusters is displayed since no less cluster can be obtained with this method. We observe that trajectories are well classified according to their path. In Fig. 11, clusters found using HCA WARD seem to be consistent. The cluster size with AP method is 36. This is a large number according to the Within-Like criterion computed with HCA. In fact, the Within-Like criterion does not decrease much between 15 and 36. However, we can see that the number of clusters found with AP are still consistent.

A cluster computed with the HCA WARD method based on a matrix distance computed with SSPD gives the best result. The Between-Like and Within-Like criteria show that this method is best used to regroup cluster around exemplar. We obtain a partition of the trajectories subset, such as each cluster represents a path taken by the drivers. We obtain a partition of traffic based on the taxi drivers' behavior leaving the Caltrain

VII. CONCLUSION

Clustering of non Euclidean objects deeply relies on the choice of a proper distance. For trajectories analysis, we presented different distances focusing on different features of such objects. To cope with their different weaknesses we propose a new distance, the Symmetrized Segment-Path Distance. This distance is time insensitive, and compares the shape and the physical distance between two trajectory objects. It does not require any additional parameters nor mapping trajectories in a different space. Hence, It can be applied on any set of trajectories, regardless of the area they come from. It enables us to obtain good clustering using either hierarchical clustering and affinity propagation methods. In this way, the clusters obtained are homogeneous with regard to shape and seem to properly capture the behaviours of the drivers. We have thus obtained a partition of the network based on the drivers' usage that can still be interpreted as vehicle trajectories. This partition can be used to solved different problem. Many applications which recommend places to visit, or which target advertising based on our destination need to forecast the final destination of drivers or predict the travel time of driver trips. Cities which wish to organize trip distribution of a city, also need to know the behaviours of the cars drivers. Some of these problems will be tackled in a following work, based on the partition obtained with our method to cluster trajectories.

Fig. 1 :

 1 Fig. 1: Supremum of P oint -to -Segment distance from point of segment s 1 1 to segment s 2 1 .

Fig. 2 :

 2 Fig. 2: Grid representation of a segment

4 )

 4 Pros and Cons:• Frechet and Hausdorff distances are both metrics, meaning they satisfy triangular inequality. With clustering algorithms like dbscan or K-medoid this is a necessary property of the distance used if we want the clustering algorithm to be efficient. They have been widely used in many domains where shape comparison is needed. But they can fail to compare trajectories as a whole. Indeed both Fréchet and Hausdorff distance return a maximum distance between two objects at given points within the two objects. As we can see in Fig.3, despite the fact that the trajectories T 1 and T 2 are well separated at the maximum value of x, they are clearly more similar to each other than to T 3 . But with a Hausdorff calculated distance, there are no strong differences between D Haussdorf (T 1 , T 2 ) = 3.26, D Haussdorf (T 1 , T 3 ) = 3.02 and D Haussdorf (T 2 , T 3 ) = 3.5. With Frechet, D F rechet (T 1 , T 2 ) = 6 is even bigger than both D F rechet (T 1 , T 3 ) = 4.19 and D F rechet (T 2 , T 3 ) = 4.17.

Fig. 3 :

 3 Fig. 3: Frechet And Hausdorff Computation between three trajectories

Fig. 4 : 1 Proposition 1 .Fig. 5 :

 4115 Fig. 4: Distance from point p 2 1 to trajectory T 1 Proposition 1. If the points that compose the trajectory T 1 lie in the set of points p 2 pl that compose the piece wise linear representation, T 2 pl , of trajectory T 2 , then D SP D (T 1 , T 2 ) = 0. Proof. If the points that compose the trajectory T 1 lie in the set p 2 pl that compose the piece wise linear representation, T 2 pl , of trajectory T 2 , all points of T 1 lie within one of the segments, that compose T 2 pl . By definition D ps (p 1 i1 , s 2 i2 ) = 0

Fig. 6 :

 6 Fig. 6: Trajectories subset

Fig. 8 :Fig. 9 :

 89 Fig. 8: Evolution of the Between-Like (a) and Within-Like (b) criteria depending on cluster size for all distances using the HCA-WARD method

Fig. 10 :

 10 Fig. 10: Clustering results with SSPD distance

Fig. 11 :

 11 Fig. 11: The isolated clusters

TABLE I

 I

		: Notation
	T	The set of trajectories
	T i	The i th trajectory of set T
	T i pl	The piece wise linear representation of T i
	n i	Length of trajectory T i
	n i pl	Length of the T i pl
	p i k	The k th location of T i
	p i pl	The set of continuous points that compose T i pl
	s i k	The line segment between p i j and p i k+1
	t i k	The time index of location p i

k p k p l 2

The Euclidean distance between p k and p l

TABLE II :

 II Metric Definition

	Property	Metric Name	d i s s i m i l a r i t y	s y m m e t r i c	m e t r i c
	Non-Negativity

TABLE III :

 III Re-Indexing based distance definition

	Cost function	Distance
	δ

TABLE IV :

 IV Re-Indexing based distance properties

	Name	Metric Types	Computation
			Cost
	DTW	symmetric	O(n 2 )

TABLE V :

 V Shape based distance properties

	Name	Metric Types	Computation
			Cost
	Hausdorff	metric	O(n 2 )
	Frechet	metric	O(n 2 log(n 2 ))
	discrete Fréchet	symmetric	O(n 2 )
	OWD	symmetric	O(n 2 log(n))
	OW D grid	symmetric	O(mn)

TABLE VI :

 VI Computation Time in seconds

	Distance	Computation time
	Fréchet	268.32
	Discrete Fréchet	0.58
	Hausdorff	2.47
	DTW	0.66
	LCSS	0.60
	SSPD	2.46
	OWD Grid 5	1.88
	OWD Grid 6	7.44
	OWD Grid 7	52.96

ACKNOWLEDGMENTS

We thank the referees for careful reading and numerous suggestions. They led us to an improvement of the work We also thank Frederic and Susan Bejina for helping improve the clarity of the paper.