Luc Brun
email: luc.brun@ensicaen.fr

Alessia Saggese

Mario Vento
email: mvento@unisa.it

Learning and Classification of car trajectories in road video by string kernels

Keywords: Abnormal Trajectories Recognition, Abnormal Behaviors Recognition, Clustering

An abnormal behavior of a moving vehicule or a moving person is characterized by an unusual or not expected trajectory. The definition of exptected trajectories refers to supervised learning where an human operator should define expected behaviors. Conversely, definition of usual trajectories, requires to learn automatically the dynamic of a scene in order to extract its typical trajectories. We propose, in this paper, a method able to identify abnormal behaviors based on a new unsupervised learning algorithm. The original contributions of the paper lies in the following aspects: first, the evaluation of similarities between trajectories is based on string kernels. Such kernels allow us to define a kernel-based clustering algorithm in order to obtain groups of similar trajectories. Finally, identification of abnormal trajectories is performed according to the typical trajectories characterized during the clustering step. The method has been evaluated on a real dataset and comparisons with other state-of-the-arts methods confirm its efficiency.

INTRODUCTION

In the last decades we have witnessed a growing need for security in many public environments, which has lead to a proliferation in the number of control systems and, consequently, in the presence of acquisition peripherals. In particular, cameras represent a suitable solution for their relative low cost of maintenance, the possibility of installing them virtually everywhere and, finally, the ability to analyze complex events.

For these reasons, a deep analysis has been recently conducted in order to realize control systems able to automatically generate alarms. Most of researches recently conducted in the field of behavior analysis has focused on the recognition of simple activities (i.e. running, waving, jumping) in high resolution videos, by exploiting the details of human body [START_REF] Aggarwal | Human activity analysis: A review[END_REF]. The main problem in such an approach lies in the fact that in a lot of real applications detailed information related, for instance, to the pose or to the clothing colors of people are not available.

As a matter of fact, in these situations, objects are in a far-field or video has a low-resolution: the only information that a video analytic system is reliably able to extract is a noisy trajectory. This consideration has recently drawn the scientific community to store (d'Acierno et al., 2012a) (d'Acierno et al., 2012b) and analyze moving objects' trajectories in order to understand their behaviors, identifying abnormal ones [START_REF] Piciarelli | Trajectory-based anomalous event detection[END_REF] [START_REF] Acampora | Combining neural networks and fuzzy systems for human behavior understanding[END_REF]. In fact, in a lot of real contexts trajectories provide the control system with enough elements to detect an anomalous behavior inside a scene: think, as an example, to a person who moves in the opposite direction of a crowd or who follows a path that the system has never seen.

The architecture of a system for behavior understanding is usually based on the following steps: learning phase and operating phase. The learning phase aims at defining rules or at extracting prototypes of normal trajectories. The definition of rules is strongly dependent on the environment and, at the same time, on the knowledge that the human operator has about the possible (ab)normal behaviors [START_REF] Dee | Detecting inexplicable behaviour[END_REF]. On the other hand, the extraction of prototypes for (ab)normal trajectories can be performed by following one of this two models: supervised and unsupervised. Techniques trained in supervised mode [START_REF] Zhou | Detecting anomaly in videos from trajectory similarity analysis[END_REF] assume the availability of a training data set with labeled instances of normal as well as abnormal trajectories. However, such an approach has a significant drawback: abnormal instances are usu-ally far fewer compared to normal ones in the training set, so implying that the prototypes extracted for abnormal trajectories are not accurate and representative. Furthermore, it is impossible to predict all abnormal behaviors inside a complex scenario.

Techniques operating in unsupervised mode [START_REF] Morris | Trajectory learning for activity understanding: Unsupervised, multilevel, and long-term adaptive approach[END_REF] do not require labeled data since they make the implicit assumption that normal instances are far more frequent than abnormal ones. An unsupervised learning phase makes the control system context-independent and can be easily applied in different real environments, since it does not use human knowledge. This is a very important and not negligible feature, since it allows the system to autonomously understand dynamics within a scene.

For all these reasons, we propose an unsupervised approach, where an abnormal trajectory refers to something that the control system has never (or rarely) seen. However, a system that raises an alarm for each trajectory which has not been seen before risks to generate too many false alarms: the system needs to identify a normal trajectory as one enough similar to one or more trajectories that the system already knows. For this purpose, we propose a learning phase based on the following steps, as depicted in Figure 1a:

• Trajectory extraction: the tracking algorithm detailed in [START_REF] Di Lascio | Tracking interacting objects in complex situations by using contextual reasoning[END_REF] is applied in order to extract moving objects' trajectories from a video for a long time period.

• Trajectories representation: the scene is partitioned into zones according to the distribution of trajectories; starting from this, each trajectory is represented as a sequence of symbols, according to the zones crossed in the scene.

• Trajectories similarity: similarity between two trajectories is evaluated by using a kernel-based method. The main advantage in this choice lies in the fact that we may combine these kernels with a large class of clustering and machine learning algorithms, which can be expressed using only scalar product between input data.

• Clustering: Given the kernel, a novel clustering algorithm is applied in order to extract clusters of trajectories inside the scene. Each cluster encodes a type of normal trajectories, dynamically extracted from the scene.

Once extracted the prototypes of normal trajectories, the control system can start the operating phase, depicted in Figure 1b: for each detected abnormal trajectory, it raises an alarm. In particular we propose to subdivide the operating phase in the following steps: • Trajectory extraction: the trajectory is extracted from a video by using the tracking algorithm detailed in [START_REF] Di Lascio | Tracking interacting objects in complex situations by using contextual reasoning[END_REF].

• Trajectory representation: the extracted trajectory is represented as a sequence of symbols.

• Classification: the trajectory is compared with the prototypes of each cluster and a similarity value is obtained for each comparison.

• Decision: the computed similarity values are processed; if such similarities are sufficiently high the trajectory is considered normal (✓), otherwise it is considered abnormal (✗). In this way, the proposed system is able to identify both rare and atypical trajectories: the former refer to something that does not appear in the training set (or only rarely appears); the latter consider all those trajectories differing in a slightly but significant way from a group of normal trajectories.

In this paper we focus on the learning phase: Subsection 2.1 shows the algorithm used to adaptively partition the scene into zones, while trajectories representation is presented in Subsection 2.2. Some details about the metric used to evaluate the similarity are provided in Subsection 2.3, while Section 3 provides a description of the proposed clustering algo-rithm. Furthermore, Section 4 shows the approach used to verify if a novel trajectory is normal or abnormal. Experimental results, which confirm the efficiency of the proposed method, are finally presented in Section 5.

THE LEARNING PHASE

A trajectory t can be seen as a sequence of k spatio-temporal points p i = [p i

x , p i y , p i t ]: t =< p 1 , p 2 , ..., p k >. This representation has two main drawbacks: first, a trajectory results in a very large amount of data to be managed; second, row data are more sensible to noise and tracking errors, and thus a filtering of each trajectory is needed before use. Furthermore, if a system considers the similarity between row data, it can introduce non relevant differences between trajectories. For example, many trajectories on a garden path may be considered as similar independently of the exact position of people on the path.

For this reason, a common representation of a trajectory consists in a reduced sequence of symbols, namely a string, aiming to preserve only the discriminant information and to reduce the space required to store trajectories.

The discriminant information to be preserved is strongly influenced by the aim of the system: as a matter of fact, in order to verify, for instance, if a person is moving in the opposite direction of a crowd or if a vehicle is driving on the emergence line on the highway, the most discriminant feature is the sequence of zones crossed by the moving object. Such scenarios can be labeled as constrained: the moving objects are expected to follow given paths within the scene.

From these considerations, we need to partition the scene into a set of zones, hence associating a single symbol to a sequence of points and eliminating non discriminant information. The criterion adopted to subdivide the scene certainly influences the performance of the entire system. As a matter of fact, on one hand it strongly affects the time needed for the computation of similarity between trajectories; on the other hand, it could decrease the reliability of the system if the chosen number of zones is not sufficiently representative of the scene. The simplest way could be to divide the scene using a fixed-size uniform grid, as in [START_REF] Papadopoulos | Trajectory retrieval with latent semantic analysis[END_REF]. The main drawback of an uniform grid, however, is that each zone has an uneven statistics, causing only a suboptimal statistical segmentation of trajectories. Furthermore, it is evident that the distribution of trajectories in the scene highlights region of interests, in which the major parts of trajectories lie and for which we need an higher level of detail.

In order to overcome these limitations, we propose an adaptive method aimed at minimizing the mean error made when assimilating a trajectory to its zone (Section 2.1). As a consequence of this partitioning criterion, areas in the scene in which most of trajectories lie are represented with an higher number of zones.

Scene Partitioning

Initially, the scene is represented by a single zone Z 1 . The scene partitioning algorithm aims at dividing Z 1 into a fixed number N of zones. The main idea behind our algorithm is to exploit the distribution of the training set by taking into account the density, as in the clustering algorithm proposed in [START_REF] Brun | Digital Color Imaging Handbook, chapter 9 : Color quantization[END_REF]. Each zone (Z i ) i∈{1,...,N} is represented by using its statistical properties (mean, major axis and covariance matrix); the proposed method works by recursively splitting a selected zone by a set of planes (cutting planes) at a chosen position (cutting position). In the following, more details about this algorithm will be provided since an enhanced kernelized version of it will be introduced in Section 3.

Let p be a generic point in the scene. We define f(p) as the number of trajectories passing through p in the image.

Using function f (.) and our zone's representation, we define the statistical property of a zone Z i (cardinality |Z i |, Mean µ i and Covariance Cov i ):

|Z i | = ∑ p∈Z i f (p) µ i = ∑ p∈Z i f (p)p |Z i | Cov i = ∑ p∈Z i f (p)p • p t |Z i | -µ i • µ t i .
(1)

A splitting strategy is based on the definition of the following steps:

• The selection of the next cluster to be split; • The selection of the cutting axis (i.e. the direction normal to each cutting plane); • The selection of the cutting position (i.e. the location of the cutting plane along the cutting axis). Since the set of all possible partitions into N zones is too large for an exhaustive enumeration, an heuristic needs to be applied. In particular, we decided to use the following heuristics, detailed in [START_REF] Braquelaire | Comparison and optimization of methods of color image quantization[END_REF].

Splitting Strategy: Each zone Z i is recursively split into two sub-zones until the N final zones are obtained. This bipartitioning strategy generates a treestructured vector quantization, as shown in Figure 2. Each leaf of the tree represents a zone of the scene. Cluster Selection: At each iteration, a single leaf of the tree-structured vector quantization is selected and then split into two zones. Therefore, the choice of the zone to be split plays a crucial role. Our algorithm attempts to minimize the total squared error T SE induced by the partition P = {Z 1 , ..., Z N }:

T SE(P) = N ∑ i=1 SE(Z i ), (2)
where SE(Z i ) is the squared error of one zone Z i computed as follows:

SE(Z i ) = ∑ p∈Z i ||µ i -p|| 2 . ( 3
)
In order to minimize T SE(P), the algorithm splits the zone Z with the maximum squared error, since its contribution to the T SE(P) is the largest. It has been shown in [START_REF] Wan | An algorithm for multidimensional data clustering[END_REF]) that this strategy corresponds to a good compromise between the computational time and the final quantization error.

Cutting Axis: Given a zone Z i to be split, we need to determine the location of the cutting plane. As in [START_REF] Braquelaire | Comparison and optimization of methods of color image quantization[END_REF], we decide to split along the axis with the greatest variance, namely the major axis.

Cutting Position: Once chosen the cutting axis, we need to select the optimal cutting position t ⋆ to subdivide the zone

Z i into two sub-zone Z 1 i ⋆ and Z 2 i ⋆ .
In particular, we choose the value able to maximize the decrease of the total squared error induced by the split:

SE(Z i ) -SE(Z 1 i ) + SE(Z 2 i )] (4)
Let m and M be respectively the minimum and the maximum projections of Z i on the cutting axis. It can be proved [START_REF] Braquelaire | Comparison and optimization of methods of color image quantization[END_REF]) that the maximization of Equation 4 can be reached by computing the optimal cutting position t ⋆ as follows:

t ⋆ = arg max t∈[m,M] δ(t) 1 -δ(t) ||µ i -µ 1 i || 2 , ( 5
)
where µ i and µ 1 i denote respectively the mean of Z i and Z 1 i ⋆ and δ(t) is defined as

δ(t) = |Z 1 i | |Z i | .
It is worth noting that if δ(t) = 1/2, the zone is divided into two sub-zones with the same cardinality.

An example of the output of the proposed scene partitioning method is provided in Figure 3: starting from the trajectory distribution in Figures 3a and3b, we show the partition of the scene by using different values of N = {20, 50}.

Trajectory representation

Once partitioned the scene into zones, a trajectory can be segmented into l segments: t = {< s 1 , ..., s l > }, where the j-th segment is defined as the sequence of points lying in the same zone Z k :

s j = {p i ∈< p a , ...p b > |p i ∈ Z k }.
Finally, the operator α(•) allows us to map the j-th segment into a symbol of our alphabet, each symbol identifying the passing through a zone. It means that the trajectory t can be now defined as t = {< α(s 1 ), ..., α(s l ) >}.

Once obtained the information about the zones crossed in the scene, additional features are extracted in order to improve the reliability of the proposed system. For each segment, information about the speed v and the shape s of the trajectories are taken into account by means of the operator θ(•). In particular, we use the Bernstein Polynomial Approximation to model each trajectory into a zone:

s = c s ∑ i=0 = a s i • t i v = c v ∑ i=0 = a v i • t i . (6)
In our experiment, c s is bounded by 3 and c v is bounded by 2. The operator θ(•) is then computed as the vector composed by the a i values:

θ(•) = [a s 1 , ..., a s c s , a v 1 , ..., a v c v ].
Thanks to this representation, each trajectory can be seen as t = {< α(s 1 ), ..., α(s l ) >, < θ(s 1 ), ..., θ(s l ) >}.

Although the operator α(•) gives to our method the most important contribution, the shape of the trajectory is an important feature, since it contributes to distinguish trajectories lying in the same zones but with very different shapes. We can think, for instance, to a person which moves with a very irregular shape but following a common path.

Trajectories similarity

The complexity and the different typology of information to take into account to represent a trajectory result in a complex strategy to verify the similarity between trajectories. In fact, we need to manage a string for the position and a sequence of numerical values for the speed and the shape, respectively obtained by means of the α(•) and the θ(•) operators.

In the last years, a lot of different methods based on dynamic programming have been proposed in order to evaluate the similarity between two sequences. These algorithms are based on similarity criteria such as the Dynamic-Time-Warping (DTW) score [START_REF] Shimodaira | Dynamic time-alignment kernel in support vector machine[END_REF], the Smith Waterman algorithm [START_REF] Saigo | Protein homology detection using string alignment kernels[END_REF] and the edit-distance [START_REF] Neuhaus | Edit distance-based kernel functions for structural pattern classification[END_REF]. However, the main problem lies in the fact that, although these methods are able to compute a similarity value, they do not define a metric. In order to solve these problems, we propose a novel similarity metric based on kernels: the main advantage is that the problem can then be formulated in an implicit vector space on which statistical methods for pattern analysis can be applied.

Furthermore, thanks to this choice, is it possible to evaluate the similarity between sequences of symbol with different length, so avoiding to force the representation of trajectories to a vector of features with a fixed dimension, as in [START_REF] Piciarelli | Trajectory-based anomalous event detection[END_REF].

In particular, we construct our kernel starting from the Fast Global Alignment Kernel (FGAK) proposed in [START_REF] Cuturi | Fast global alignment kernels[END_REF]. The main idea of all global alignment kernels is to measure the similarity between two sequences by summing up scores obtained from local alignments with gaps of the sequences.

An alignment between two sequences x = {x 1 , ..., x n } and y = {y 1 , ..., y m } of length n and m respectively is a pair of increasing integral vectors

(π 1 , π 2 ) of length p < n + m, such that 1 = π 1 (1) ≤ ... ≤ π 1 (p) = n and 1 = π 2 (1) ≤ ... ≤ π 2 (p) = m,
with unary increments and no simultaneous repetitions. Let A(n, m) be the set of all the possible align-ments between the two time series of lengths n and m.

The global alignment kernel (GAK) is defined as:

k GA (x, y) = ∑ π∈A(n,m) |π| ∏ i=1 k(x π 1 (i) , y π 2 (i) ). (7)
It can be shown [START_REF] Cuturi | Fast global alignment kernels[END_REF]) that k GA is a positive definite kernel if k and k/(1 + k) are positive definitive kernels. Furthermore, the GAK avoids the diagonal dominance of the Gram matrix. Diagonal dominance is an undesirable property, since it implies that all the points in a training set are nearly orthogonal to each other in the corresponding feature space.

Starting from the representation of our trajectories, we need to define a kernel which is able to properly combine all the different features related to a trajectory. In particular, we defined the following kernels.

Toeplitz Kernel: In order to speed up the computation of the kernel, we use the triangular kernel for integers, also known as Toeplitz kernel:

w(i, j) = 1 - |i -j| T , ( 8
)
where T is the order of the kernel. The main advantage in the use of the triangular kernel is that it allows to only consider a smaller subset of alignments. Dirac Kernel: In order to evaluate the similarity between two strings α(x) and α(y) encoding the sequences of zones respectively traversed by trajectories x and y, we use a dirac kernel δ(α(x i ), α(y i )), defined as:

δ(α(x i ), α(y i )) = 0 if α(x i ) = α(y i ) 1 if α(x i ) = α(y i ) (9)
The Dirac Kernel is combined with the Toeplitz Kernel so obtaining:

k Z (x i , y j , i, j) = w(i, j) • δ(α(x i ), α(y j )). ( 10
)
Weighted Dirac Kernel: The main lack of this similarity evaluation lies in the fact that the proximity of two zones is not considered. In order to overcome this limitation by taking into account adjacency relationships between zones, a weighted version of the dirac kernel is also exploited:

k W Z (x i , y j , i, j) = w(i, j) • δ w (α(x i ), α(y j )). (11)
Zones are mapped into a non-oriented weighted graph G = {V, E, w}, whose vertices V = {V 1 , ...,V N } identify zones and whose edges E = {E 1 , ..., E L } identify proximity of two zones. Each edge is associated to a weight e v 1 ,v 2 , identifying the number of pixels separating two zones. An example is shown in Figure 4:

δ w (α(x i ), α(y i )) =          0 if α(x i ) = α(y i ) and e α(x i ),α(y i ) / ∈ E e I α(x i ),α(y i ) if e α(x i ),α(y i ) ∈ E 1 if α(x i ) = α(y i ) (12
) where e I α(x i ),α(y i ) is a normalized version of e α(x i ),α(y i ) , obtained by dividing e α(x i ),α(y i ) by the length of the longest zone's border.

Speed and Shape Kernel: The evaluation of the similarity related to the velocity and to the shape is based on the following kernel:

k SS (θ(x i ), θ(y i )) = e -φ σ (θ(x i ),θ(y i )) , (13)
where

φ σ (θ(x i ), θ(y i )) = 1 2σ 2 ||θ(x i ) -θ(y i )|| 2 + log 2 -e - |θ(x i )-θ(y i )|| 2 | 2σ 2 . ( 14
)
This last kernel is used instead of the Gaussian one in order to guarantee the p.d. of k GA [START_REF] Cuturi | Fast global alignment kernels[END_REF]. The combination of these two last kernels is defined as:

k (W )ZSS (x i , y j , i, j) = k (W )Z (α(x i ), α(y j )) • k SS (θ(x i ), θ(y i )). ( 15
)
Starting from Equation 7, products of any of the 4 kernels (k Z , k W Z , k ZSS and k W ZSS ) can be considered to obtain the final kernel k GA . Finally, a normalization of the kernel is performed in order to normalize kernel's values in the interval [0, 1]. Therefore, the final normalized kernel k N GA is:

k N GA (x i , y j , i, j) =
k GA (x i , y j , i, j) k GA (x i , x i , i, i) * k GA (y j , y j , j, j) .

(16)

CLUSTERING

From a general point of view, the goal of a clustering algorithm is to find a fixed number N C of groups that are both homogeneous and well separated, that is, trajectories within the same group should be similar and entities in different groups dissimilar. In our context, we aim at exploiting a clustering algorithm in order to obtain a set of prototypes of normal trajectories.

In the last decades, a lot of graph-based clustering algorithms [START_REF] Foggia | Assessing the performance of a graph-based clustering algorithm[END_REF] [START_REF] Schaeffer | Graph clustering[END_REF] [START_REF] Foggia | A graph-based algorithm for cluster detection[END_REF], like Spectral Clustering or Cut-Clustering, have been exploited. Although these techniques seem to provide good results, they do not allow to readily verify if a novel trajectory belongs to a cluster, that is our main objective.

In order to avoid this restriction, k-means approach and its derivative methods are most frequently used. In particular, the Kernel k-mean [START_REF] Dhillon | Kernel kmeans: spectral clustering and normalized cuts[END_REF]) is a generalization of the standard k-means algorithm: the input data are mapped into a higher dimensional feature space through a non-linear transformation and then k-means is applied in feature space. In this way, this algorithm allows to separate the non linearly separable clusters. However, the main problem of such an approach consists in the initialization, which strongly influences the performance of this method since the algorithm converges to the local minimum closest to the initial condition.

In [START_REF] Tzortzis | The global kernel kmeans algorithm for clustering in feature space[END_REF] an improved version of the basic Kernel k-means, the Global Kernel k-means, has been proposed. The main idea is that a near-optimal solution with k clusters can be obtained by starting with a near-optimal solution with k -1 clusters and initializing the kth cluster appropriately based on a local search. During the local search, N initializations are tried, where N is the size of the data set. The k -1 clusters are always initialized to the k -1-clustering problem solution, while the kth cluster for each initialization includes a single point of the data set. The solution with the lowest clustering error is kept as the solution with k clusters. Since the optimal solution for the 1-clustering problem is known, the above procedure can be applied iteratively to find a near-optimal solution to the M-clustering problem. It is clear that the main drawback of the Global Kernel k-means is the high computational cost, as experimental results conducted will show in Section 5.

In order to overcome these limitations, we introduce a novel and efficient kernelized clustering algorithm. The proposed clustering algorithm is based on the splitting methods presented in Section 2.1: the cluster with the maximum squared error is selected and then split into two different clusters along the major axis. However, the main novelty refers to the kernelization of the considered algorithm.

Thanks to the chosen heuristics, the partitioning of the space into N C clusters is performed by a sequence of N C -1 iterations. It is an important and not negligible feature, since a lot of recently proposed clustering algorithms [START_REF] Tzortzis | The global kernel kmeans algorithm for clustering in feature space[END_REF] are very expensive from a computational point of view, as we will show in Section 5.

Given the cluster C containing all the trajectories belonging to the training set, let us consider a generic cluster C t ⊂ C. This cluster is encoded by the following vector of R |C| :

Λ i C t = 1 if i ∈ C t 0 otherwise (17)
Finally, K is the Gramm Matrix of the training set, defined by K i j = (k N GA (s i , s j )) i j . As mentioned in Section 2.1, our clustering algorithm builds a sequence of partitions P 1 , . . . , P n of C, with P 1 = {C} and P n enconding a partition of C into n clusters. The following heuristics have been selected for each step:

Cluster Selection: For each iteration k, the cluster C t of P k with the maximum squared error SE(C t ) is selected. In particular, we decided to cut the cluster along the major axis, obtained by means of a Kernel PCA [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF]. The Gramm Matrix K is first diagonalized and the centered matrix

SE(C t ) = ∑ s∈C t ||ψ s -µ t || 2 = |C t | - 1 |C t | Λ t C t KΛ C t , (
K ′ = (K -1 C t • K -K • 1 C t + 1 C t • K • 1 C t ) is obtained, with 1 C t being a |C t | by |C t | matrix for which each element takes value 1/|C t |.
The first eigenvector α, satisfying: λα = K ′ α is then computed. For each trajectory s, the projection of ψ s on the major axis α is obtained by:

< α, ψ s >= |C t | ∑ i=1 α i k(s i , s). ( 19
)
Thanks to the computed projections, trajectories belonging to C t are ordered along the major axis in order to obtain the optimum cutting position.

Cutting Position: Given the selected cluster C t and its cutting axis, the cutting position t ⋆ is computed in the range [m, M], being m and M respectively the minimal and the maximal projection of C t on the major axis by means of Equation 5.

It can be shown A.2 that ||µµ t || 2 can be computed as follows: We can note that this operation must be performed for each point in the range [m, M]. Since it requires multiple matrix multiplications, it results in a high computational cost. Let p denotes the next trajectory to add to C t in order to obtain C t+1 (C t+1 = C t ∪ {p}). It can be shown A.2 that ||µ-µ t+1 || 2 can be efficiently updated from ||µµ t || 2 . In particular, the first term 1 |C| 2 e t Ke is constant since it does not depend on the partition induced by t. Let be Λ C t+1 = Λ C t + δ p , being δ p the vector of zeros containing a single 1 at position p. The second term e t KΛ C t+1 of equation 20 may be defined iteratively as follows:

||µ-µ t || 2 = 1 |C| 2 e t Ke- 2 |C t ||C| e t KΛ C t + 1 |C t | 2 Λ t C t KΛ C t , (
e t KΛ C t+1 = e t KΛ C t + n ∑ i=1 k(i, p). ( 21
)
Note that the second term of equation 21 may be precomputed. Equation 21 is thus evaluated in constant time. Finally, the last term Λ t C t+1 KΛ C t+1 becomes:

Λ t C t+1 KΛ C t+1 = Λ t C t KΛ C t + 2 ∑ i∈C t k(i, p) + k(p, p).
(22) Using equations 21 and 22, we significantly reduce the computational cost of our algorithm: as a matter of fact, the evaluation of µµ t+1 2 only requires to compute, for each iteration, values ∑ i∈C t k(i, p) and k(p, p) this last term being equal to 1 since we use a normalized kernel.

Stop Condition: In our context, the clustering algorithm is used to initialize the system during its unsupervised learning phase. It means that the number of clusters can not be fixed a priori, since the system has not any knowledge about the environment. For this reason, we choose to use as stop condition a lower bound on the mean squared error (MSE) made when assimilating one trajectory to its cluster:

MSE(C t ) = SE(C t ) |C t | . ( 23
)
In this way, the system does not need knowledge of the human operator about the environment, but is able to determine the optimum number of clusters starting from the distribution of trajectories.

OPERATING PHASE

The operating phase aims at identifying abnormal behaviors according to the set of typical trajectories determined during the learning phase (Section 3). In particular, our algorithm evaluates the distance between a trajectory t s and all cluster's centers C 1 , ...C N C obtained during the learning phase. The cluster with the closest mean from t s is selected as the potential typical trajectory followed by t s . An additional test should then be performed in order to determine if t s belongs to this closest cluster. According to this last test t s is classified as normal (it belongs to one cluster encoding typical trajectories) or an alert is raised and t s is classified as an abnormal behavior.

Classification: Les s denote the string associated to t s and ψ s the projection of s into the Hilbert space encoded by one of our kernel. The squared distance between ψ s and the mean µ t of a cluster C t is defined by:

d 2 t (µ t , Ψ s ) =< µ t , µ t > + < Ψ s , Ψ s > -2 < µ t , Ψ s > = 1 + 1 -2 < µ t , Ψ s > = 2(1-< µ t , Ψ s >) = 2 1 - 1 |C t | ∑ s i ∈C t k(s, s i ) . (24)
Decision: Let C * t denote the cluster with the closest center (µ * t ) determined according to equation 24. Since our clustering algorithm always split clusters according to their axis of greatest variance, we consider that the covariance matrix of each cluster is approximately diagonal. In this case, a threshold on the Gaussian probablity that string t s belongs to C * t is approximated by comparing the squared distance d 2 (µ * t , Ψ s ) with a multiple of the squared error of C * t :

d 2 (µ * t , Ψ s ) ≤ α * MSE(C * t ).
(25) Conversely to the parameter ν of one class SVM [START_REF] Cortes | Support-vector networks[END_REF], a high value of α provides a better generalization but may increases the number of false positive in the test determining the classification to C * t .

EXPERIMENTAL RESULTS

The proposed method has been validated on the MIT trajectories dataset [START_REF] Wang | Trajectory analysis and semantic region modeling using nonparametric hierarchical bayesian models[END_REF], a standard and freely available dataset composed by 40.453 trajectories obtained from a parking lot scene within five days. Starting from the entire dataset D, a subset D * of trajectories belonging to vehicles (10.335) has been manually extracted by an expert (see Figure 3a) and the proposed system has been evaluated. The experiments have been conducted on a MacBook Pro equipped with Intel Core 2 Duo running at 2.4 GHz.

The experimentation has been conducted into three different steps: the first aims to validate the proposed kernels, highlighting advantages and disadvantages of each by means of a visual interpretation. The second is a quantitative and qualitative experimentation conducted over the clustering algorithm, aimed at confirming the efficiency of the proposed method if compared with other state of the art methods. Finally, a misclassification matrix is presented, in order to confirm the efficiency of the proposed classification method.

Experimentation over kernels: As for the first step, it is interesting to note that kernels induce metrics. This means that a quantitative comparison of different kernels is not possible outside a regression or classification task with a ground truth which is not available for this data set. We thus only perform a qualitative evaluation based on visual comparisons. Such an experiment has been performed over the dataset D * , aiming at extracting the M most similar trajectories with respect to the one shown in Figure 5a (M=30).

The parameters of the kernels that we need to tune are respectively T (for the Toeplitz Kernel) and σ (for the Speed and Shape Kernel). In particular, in [START_REF] Cuturi | Fast global alignment kernels[END_REF] is shown that the best result can be achieved by using the following parameters:

• The Bandwidth σ is set to a multiple of a simple estimate of the median distance K of different points observed in different time-series of our training set, scaled by the square root of the median length of time-series in the training set. In particular, in [START_REF] Cuturi | Fast global alignment kernels[END_REF] is suggested to try σ ∈ {0.1, 1, 10} • K, where

K = median||θ(x i ) -θ(y i )|| median(||x||) (26)
and to use higher multiples (e.g. 2,5). Being K = 60 in our dataset, we choose σ = 2 * 60 = 120. • The Triangular parameter T can be set to a reasonable multiple of the median length, e.g. 0.2 or 0.5. Being the median length equal to 5 in our dataset, we chose T = 3.

An example of our different kernels at work is shown in Figure 5. It is worth noting that the Weighted Dirac Kernel (Figure 5c) s less sensitive than the Dirac Kernel (Figure 5b) to small variations in the position of trajectories. It is a very desirable property, especially in the operating phase, since it allows to reduce the number of false positive. On the contrary, the combination of Weighted Dirac Kernel, Gaussian and Toeplitz ones (Figure 5d) increases the reliability of the evaluation by taking into account the shape of trajectories.

Experimentation over clustering: As for the second step, a qualitative evaluation of the proposed method over the dataset D * is shown in Figure 6, where some of the obtained clusters are depicted in a tree structure.

In order to have a quantitative measure, we compute the C-index [START_REF] Hubert | Quadratic assignment as a general data analysis strategy[END_REF]. It is defined as:

C = S -S min S max -S min , ( 27
)
where S is the sum of distances over all pairs of objects form the same cluster, n is the number of those pairs and S min is the sum of the n smallest distances if all pairs of objects are considered. Likewise S max is the sum of the n largest distances out of all pairs. The C-index ranges from 0 to 1 and the optimum value is 0. The above mentioned index has been used to compare the proposed method with other state of the art approaches. In particular, we considered the traditional Kernel k-mean with its two improved versions, the Global Kernel k-means [START_REF] Tzortzis | The global kernel kmeans algorithm for clustering in feature space[END_REF] and the Fast Global Kernel k-means [START_REF] Tzortzis | The global kernel kmeans algorithm for clustering in feature space[END_REF].

A comparison in terms of C-Index and computational cost is shown in Tables 1a and1b. Note that the results of the Kernel k-mean (in terms of C-index and time) is obtained by taking the minimal c-index and overall time over 200 different trials, in order to limit the dependency of the results from a particular initialization.

The value that we obtained, equal to 0.0580 if the similarity value is computed by using the Dirac Kernel, while it is equal to 0.0799 when using a Weighted Dirac Kernel instead of a simple Dirac Kernel. In both cases, the efficiency of our technique is confirmed, since the proposed method outperforms the other considered approaches. Furthermore, the computational cost is not too expensive if we consider that this algorithm does not need to work in real time, since it is only used at the start-up of the system. Experimentation over operating phase: The dataset D * has been divided into three folds and one of these has been used for the learning phase. The remaining two folds have been mixed with the remaining trajectories (D \D * ) and are used to test the system. The tests have been performed by computing the similarity between trajectories by using the Dirac Kernel and the Speed and Shape Kernel.

The confusion matrix obtained is depicted in Table 2. The results, for a fixed value of α (α = 2), show that the Weighted Dirac Kernel provides a better generalization than the Dirac Kernel, without paying in terms of false positive errors.

Furthermore, it is worth pointing out that in this case it has not been possible to compare our approach with other state-of-the-art methods since the results . . .

Figure 6: A part of the tree obtained by using the proposed clustering method. The color of each trajectory highlights its direction: it starts in green and progressively turns its color into red.

on the considered dataset are not available. As a matter of fact, a lot of recently proposed methods have been tested over owner datasets, not always available.

In any case, starting from the obtained results, which are sufficiently good for most practical applications, we can enforce the effectiveness of the method by drawing some considerations about the nature of the errors; first of all, as shown in Figure 7a, we can see a trajectory labelled as abnormal in our dataset, but classified as genuine by the system. This happens because they refer to vehicle's trajectories partially located in the grass (labelled as abnormal in the dataset) and at the same time very similar in shape and close to others that avoid the grass just for a few centimeters. The system, based for its nature on shape similarity, has no chance to give a correct answer. It is worth noting that this error cannot be recovered by any unsupervised system based on similarity evaluation; only the introduction of information about the boundaries of the areas can provide a system the possibility of discriminating among very close trajectories on the basis of possible boundary crosses. Another similar situations occurs in Figure 7b, where the vehicle tries to park, but because of car place lack, fails, consequently leaving the parking by exiting after a complete turn. In this case, the description of the trajectory follows a regular and normal path, except for a very limited stretch, reproducing the same typology of error occurring in the previous case.

Opposite kinds of error occur in Figure 8, in which some trajectories, normal with respect to their semantic and consequently labelled in this way, are wrongly classified as abnormal by the system because their short dimension, probably caused by a tracking algorithm error, do not allow to associate them with a sufficient reliability to a cluster containing normal trajectories.

In conclusions, the yet sufficiently high performance are even more acceptable at the light of the fact that a rather high percentage of errors cannot be recovered by using unsupervised methods, as their elimination would require further information possibly given by a supervised approach. However, we could think, as future work, the introduction of a mixed solution based both on clustering and boundary-constraints so as to catch the advantages of both these approaches, even at the cost of introducing a little more heavy a priory knowledge about the scene to be processed.

CONCLUSIONS

We propose a system able to identify abnormal trajectories without the human operator explicitly defining the rules. It has been achieved by introduc- ing an unsupervised method able to deduce properties of a scene from a set of trajectories. Starting from a set of normal trajectories acquired by means of a video analysis system, our method represent each trajectory by a sequence of symbols associated to relevant features of trajectories (crossed zones, shape and speed in each zone). This quantization is obtained by partitioning the scene into a fixed number of adaptive zones. Similarity between trajectories is evaluated by means of a fast alignment global kernel. Trajectories are then grouped into homogenous clusters encoding normal trajectories. The classification into (ab)normal trajectories is performed by taking advantaging on the statistical properties of the clusters. Experiments have been performed on a real dataset and the obtained results, compared with other state of the art methods, confirm the efficiency of the proposed approach. However, a deeper analysis will be conducted in order to confirm the efficiency of the entire approach if compared with other state-of-the art methods. Furthermore, some more reliable classifiers will be exploited in order to increase the overall dependability of the system.

A APPENDIX

In the following, the formal proofs of the computation of the squared error and of the cutting position of the clustering algorithm are provided.

A.1 Squared Error

Let be:

b i s =      1 -1 |C t | if i = s -1 |C t | if i ∈ C t ∧ i =
µ t = 1 |C t | ∑ i∈C t Ψ i = 1 |C t | ∑ i∈C Λ i C t Ψ i . ( 33
) b i = 1 |C| - Λ i C t |C t | ; (34) b = 1 |C| e - Λ i C t |C t | , (35)
being e a vector of one values.

It can be shown that ||µµ t || 2 can be computed as follows:

||µ -µ t || 2 = || ∑ i∈C Ψ i |C| - Λ i C t Ψ i |C t | || 2 (36) = || ∑ i∈C 1 |C| - Λ i C t |C t | Ψ i || 2 = || ∑ i∈C b i Ψ i || 2 = ∑ i, j∈C b i b j < Ψ i , Ψ j > = ∑ i, j∈C b i b j • k(s i , s j ) = b t Kb = 1 |C| e - Λ i

Figure 1 :

 1 Figure 1: Learning phase (a) and operating phase (b).

Figure 2 :

 2 Figure 2: Tree-structured vector quantization obtained recursively partitioning the scene into N = 8 zones.

Figure 3 :

 3 Figure 3: Partition of the scene starting from the training set depicted in (a) represented by the frequency map in (b). The quantization algorithm is applied with different values of N: (c) N=20, (d) N=50.

Figure 4 :

 4 Figure 4: Graph-based representation of the scene. Each zone is represented by a vertex and each border between two zones is represented by an edge. The weight of each edge is determined by the length of the corresponding border.

  20)where |C| denotes the cardinality of cluster C.

Figure 5 :

 5 Figure 5: Most similar trajectories to the one shown in (a), obtained by using Dirac Kernel (b), Weighted Dirac Kernel (c), Speed and Shape Kernel (d).

Figure 8 :

 8 Figure 7: Abnormal trajectories classified as normal.

  18) where ψ s is the projection of the string s in the Hilbert space implicitly defined by k N GA . Equation 18 may be evaluated in O(|C t | 2 ), with |C t | denoting the cardinal-

	ity of C t A.1.
	Cutting Axis: Once selected the cluster C t , we need to chose a cutting plane to obtain clusters C 1 t and C 2 t : the optimum cutting plane aims at minimizing SE(C 1 t ) + SE(C 2 t ).

Table 1 :

 1 Comparison of the proposed method with other state of the art approaches, in terms of C-index and computational cost. The similarity between trajectories is computed by using the Dirac Kernel (a) and the Weighted Dirac Kernel (b).

	Method (Dirac Kernel) C-Index Time (secs)
	Proposed Method	0.0580	947.18
		K-Means	0.3907	15.946
	Global K-Means	0.1282	3.45 • 10 4
	Fast Global K-Means	0.1877	67.882
		(a)
	Method (Weighted Dirac Kernel)	C-Index Time (secs)
	Proposed Method	0.0799	1055.98
		K-Means	0.4203	12.783
	Global K-Means	0.1376	3.45 • 10 4
	Fast Global K-Means	0.2868	75.201
		(b)
			Predicted Class
		Normal Abnormal
	GT	Normal Abnormal 15.90% 84.10%	21.40% 76.60%
		(a)
			Predicted Class
		Normal Abnormal
	GT	Normal Abnormal 14.70% 85.30%	7.10% 92.90%
		(b)
	Table 2: Misclassification Matrix obtained by using Dirac
	Kernel (a) and Weighted Dirac Kernel (b).

  ||C| e t KΛ C t + 1 |C t | 2 Λ t C t KΛ C t .Furthermore, starting from ||µµ t || 2 , ||µµ t+1 || 2 can be quickly computed. In particular, the second term e t KΛ C t+1 is:e t KΛ C t+1 = e t K(Λ C t + δ) (37) = e t KΛC t + e t Kδ p = e t KΛ C t + Finally, the last term Λ t C t+1 KΛ C t+1 becomes: Λ t C t+1 KΛ C t+1 = (Λ C t + δ p ) t K (Λ C t + δ p )

		C t |C t |	t	K	1 |C|	e -	Λ i C t |C t |
	=	1 |C| 2 e t Ke -	2 |C t n
						i=1 ∑	k(i, p).
							(38)
		= Λ t C t KΛ C t + 2Λ t C t Kδ p + δ t p Kδ p
		= Λ t C t KΛ C t + 2 ∑ i∈C t	k(i, p) + k(p, p).