# Fast maneuver recovery from aerial observation: trajectory clustering and outliers rejection

Nelson de Moura, Augustin Gervreau-Mercier, Fernando Garrido, Fawzi

Nashashibi

► To cite this version:

Nelson de Moura, Augustin Gervreau-Mercier, Fernando Garrido, Fawzi Nashashibi. Fast maneuver recovery from aerial observation: trajectory clustering and outliers rejection. 2024 IEEE Intelligent Vehicles Symposium (IV), ITS-IEEE, Jun 2024, Jeju, South Korea. hal-04497008v2

HAL Id: hal-04497008

https://inria.hal.science/hal-04497008v2

Submitted on 2 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Fast maneuver recovery from aerial observation:
trajectory clustering and outliers rejection

Nelson de Moura<fnref n="1" />, Augustin Gervreau-Mercier<fnref n="1" />,<fnref n="2" />, Fernando Garrido<fnref n="3" /> and Fawzi Nashashibi<fnref n="1" />

*Abstract*— The implementation of road user models that realistically reproduce a credible behavior in a multi-agent simulation is still an open problem. A data-driven approach is proposed here to infer behaviors that may exist in real situation to obtain different types of trajectories from a large set of observations. The data, and its classification, could then be used to train models capable to extrapolate such behavior. Cars and two different types of Vulnerable Road Users (VRU) will be considered by the trajectory clustering methods proposed: pedestrians and cyclists. The results reported here evaluate methods to extract well-defined trajectory classes from raw data without the use of map information while also separating "eccentric" or incomplete trajectories from the ones that are complete and representative in any scenario. Two environments will serve as test for the methods develop, three different intersections and one roundabout. The resulting clusters of trajectories can then be used for prediction or learning tasks or discarded if it is composed by outliers.

<figure><img src="image_2.png" /><figcaption>(a) Pedestrian on sce. 2 of inD</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(b) Vehicles on sce. 1 of inD</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(c) Vehicles on sce. 2 of rounD</figcaption></figure>

<caption>Fig. 1: Trajectory cluster for different road users</caption>

# I. INTRODUCTION

Simulation is a indispensable tool to prove the efficacy and viability of any framework or system capable to drive an Automated Vehicle (AV) before integration with a prototype. In most cases these simulations reproduce the behaviors of other road users based on real recordings, which in this case the behavior of all road users is fixed ahead of time, or it relies on hybrid approaches, mixing real information with some *a priori* hypothesis, modeling and/or knowledge about the agents being represented. The main goal of this work is to produce classifications of trajectories to feed these hybrid methods with reliable and diverse set of trajectories, for different situations and different road users while being fast to execute and reliable enough to sift through outliers at input.

Trajectory clustering has been a long research topic in the AV area. Many articles deal specially with the analysis of vehicle trajectory as a way to retrieve the possible trajectories in an environment, to study the traffic flow intersections [1], to discover possible longitudinal behaviors of vehicles [2], to execute some learning task [3] or even to examine the scenarios that might happen during driving [4]. Thus, the goal of this paper is to propose a fast and robust way to recover sets of trajectories for vehicles and vulnerable road

users (VRU), delivering sets of trajectory samples as input to all the aforementioned tasks in a simple manner.

Approaches based on clustering with dynamic time warping (DTW) are the norm in the literature. In [4] three different threshold comparisons were made using the DTW distance metric to establish a similarity relationship between scenarios involving multiple road users. K-means and fuzzy c-means were used for [1], [5] with longest common sub-sequence (LCSS) to cluster trajectories in intersections so to derive insights about the traffic flow in multiple lane cross-intersections. Changing from the urban to aerial traffic, [6] proposed a method to combine k-means with outlier removal based on information theory by the minimization of the holoentropy and achieving good results clustering flight trajectories. With a different motivation, [2] implemented a fast k-means clustering method only for vehicles and bypassing the outliers problem. And on a totally different scale [3] applied the same idea of trajectory clustering but on a city scale, learning an embedding to simplify the trajectory representation and then clustering the projected data to find vehicles with similar behavior.

The contribution of this paper is two-fold: First, to propose a fast clustering method for maneuver retrieval from real observations that do not need any map information and that is compatible with vehicles and VRUs. Second, to adapt this model to deal with trajectories that can be considered as outliers, separating these "eccentric"<fnref n="1" /> and/or erroneous instances without disturbing the clustering process. A preview of the results obtained can be seen in Figure 1. Given the

<footnote n="1">¹Nelson de Moura, Fawzi Nashashibi and Augustin Gervreau-Mercier are with INRIA, 75012 Paris, France, {nelson.demoura; fawzi.nashashibi}@inria.fr</footnote>

<footnote n="2">²Augustin Gervreau-Mercier is also with Université Paris Dauphine-PSL, 75775 Paris, France, augustin.gervreau@dauphine.eu</footnote>

<footnote n="3">³Fernando Garrido is with Valeo DSW team, 94000 Créteil, France, fernando.guarrido@valeo.com</footnote>

<footnote>This research has been funded by the plan "France Relance", grant agreement number ANR-21-PRRD-0005-01</footnote>

<footnote n="1">¹Definition on section II</footnote>
simplicity of the approach, the method presented here can be a valuable addition to the plethora that already exists.

All data used originate from a microscopic observation of road environments by a drone [7], [8], which retain a good amount of information about the behaviors of each road users, in comparison with [9] that make observations on a larger scope and well structured environment (4-lane road with signals), and of [10] that focus more on the interaction size with short observation periods (4s) and 10Hz of acquisition frequency (in comparison with 25Hz of the data used here).

## II. SEPARATING TRAJECTORIES OF INTEREST

Differently from the vehicle trajectories studied in [2], VRU trajectories are less constraint by its environment and are also more prone to acquisition error as well (shadows, changes in direction, multiple users close by). Also, when a scenario for observation is defined some of the less intuitive trajectories become superfluous, considering the interest of keeping only those which represent behaviors that can be transposed in other scenarios. Take, for example, the trajectories displayed in Figure 2 (white cross represents the beginning of the trajectory and black cross the end): these three different sets may represent a real-life situation, like getting out of a store and entering in a car but they are not of interest since they are scenario-specific. These types of trajectories will be qualified as eccentric from now on.

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 2: One eccentric trajectory and a target cluster (background image from [7])</figcaption></figure>

The main focus is to sift through entire datasets and isolate the eccentric trajectories (like the ones in Figure 2a) and erroneous ones (like vehicles starting their trajectory in the middle of an intersection) in specific clusters and trajectories of interest Figure 2b in their own clusters. The final result can then be visually inspected to discard some and retain others. Both figures were produced using the data available in the InD Dataset [7].

Given the difficulties that methods of the similar inspiration of k-means have with outliers, other clustering algorithms were considered to deal with pedestrian and cyclist trajectories. All of them using a pre-calculated dissimilarity matrix (Equation (1)) where each element is the result of the Dynamic Time Warping (DTW, subsection II-A) distance measure of two trajectories.

$$D_{DTW} = \begin{bmatrix} d_{0,0}^{DTW} = 0 & \cdots & d_{0,n}^{DTW} \\ d_{1,0}^{DTW} & \cdots & d_{1,n}^{DTW} \\ \vdots & \ddots & \vdots \\ d_{n,0}^{DTW} & \cdots & d_{n,n}^{DTW} = 0 \end{bmatrix} \quad (1)$$

### A. Dynamic time warp (DTW)

DTW was first introduced in the speech processing domain as a way to compare two time series that have different phases. Even though the trajectories studied here were sampled at the same frequency, they may have different lengths and although corresponding to the same maneuver in an intersection. Consider two discrete time series, represented by (2) and (3), with different sizes $n$ and $m$, where $K = \{k_0, k_1, \dots, k_n, \dots, k_m, \dots\}$ represents the sampled periods:

$$R[K] = r[k_0], r[k_1], \dots, r[k_n] \quad (2)$$

$$S[K] = s[k_0], s[k_1], \dots, s[k_m] \quad (3)$$

The goal of the DTW is to calculate the optimal sequence of pairs of point indexes, one from each time series. This is done by minimizing the euclidean distance between the points indicated by the index pair, from $(r[k_0], s[k_0])$ to $(r[k_n], s[k_m])$, using a certain set of increments to walk from the former to the latter. In the standard implementation (equation 5) three steps are tested: +1 on the index 1, +1 on the index 2 or +1 in both. Equation (4) defines the DTW from R and S as the calculated sum of distances, which are determined by the recursion in Equation (5), for $0 \le i \le k_n$ and $0 \le j \le k_m$.

$$DTW(R, S) = \gamma(k_n, k_m) \quad (4)$$

$$\gamma(i, j) = d(r[k_i], s[k_j]) + \min[\gamma(i-1, j), \gamma(i, j-1), \gamma(i-1, j-1)] \quad (5)$$

There are multiple DTW variants, some changing the walk used in the recursion (5) (constraint DTW [11]) or adopting restriction on the elements to be considered by Equation (5) (Sakoe-Chiba band [12]; Itakura parallelogram [13]). Usually, when the Euclidean metric is used, the centroid of a set of series can be calculated simply by summing all the elements and dividing by the number of series in the set. One can do the same with series of different lengths using the Dynamic Barycenter Averaging (DBA) [14], however the result of this algorithm usually is a non-differentiable array always with the same size of the biggest array in the set. When a time-series is necessary to represent the ensemble of a cluster, the medoid will be chosen, according to Equation (6), where $\mathcal{X}$ is the set being considered and $d$ in our case is the DTW distance.

$$x_{\text{med}} = \arg \min_{x \in \mathcal{X}} \sum_{i=0}^{N} d_{\text{DTW}}(x, x_i) \quad (6)$$
## B. Clustering methods

Three main methods were used to cluster the trajectories using the DTW distance: hierarchical clustering, partition around medoids (or k-medoids) and dissimilarity matrix clustering.

1) **Hierarchical clustering:** The hierarchical clustering used was based on a agglomerative processes, i.e., it starts with each sample being a cluster and at each step it merges the two most similar clusters, to then continue this process until the desired number of clusters is achieved [15]. The metric used to measure the similarity of two clusters, thus to decide with clusters should be merged at a given iteration, was the average linkage, Equation (7):

$$d_{C_i, C_j} = \frac{1}{N \cdot M} \sum_{x_i \in C_i}^{N} \sum_{x_j \in C_j}^{M} d(x_i, x_j) \quad (7)$$

Where $C_i$ and $C_j$ are two clusters being evaluated and $N$ and $M$ are the number of elements inside each respective cluster. The distance measure used is the DTW (from (1)).

2) **Partition around medoids (or k-medoids):** It uses the same sequence of calculation - allocation of elements in cluster then center recalculation - from the k-means algorithm but using the medoid element as the cluster center (equation 6), not a synthetic average of elements [16]. Such adaptation is common in cases where it is difficult to calculate the average of elements being clustered, as for example when these do not have the same length.

3) **Dissimilarity matrix clustering:** The algorithm was proposed to accelerate the clustering of vehicle's trajectories in [2], in comparison with the k-means algorithm. To simplify the k-means calculation, it is applied at each row of the dissimilarity matrix for the entire dataset to look for the smallest cluster center of all rows (which will be the one that has the smallest sum of distances to the element represented by the row). Then, all the elements assigned to this minimal cluster are removed from the matrix and the process is executed again, until the desired number of clusters is achieved. If there are any left, then they are assigned to the cluster in which it has the smallest distance to its medoid.

## III. METHODS EVALUATED

Independently from the clustering method used from the proposed in the previous section, some errors in classification might still appear, even if a higher number of clusters is used. To cluster using the DTW distance considers only the shape of the trajectory, which might mix trajectories that have small but important differences at its origin or terminus but that share an important part of its path. Hence, to correct this errors the initial and final points shall be used in a separate clustering procedure split the elements based on these points. Then, it will be necessary to check if any of the just-obtained sub-clusters should be fused back together, if they really are on the same maneuver, or even if they should be merged with other sub-clusters, to determine the final result. Algorithm 1 shows how the entire clustering process with this post-processing operation

works. The interval $[nk_{min}, nk_{max}]$ refers to the minimal and maximal number of clusters to be evaluated.

<table><caption>**Algorithm 1: Proposed clustering procedure**</caption><tr><td>1</td><td>for</td><td>n<sub>k</sub> &isin; [nk<sub>min</sub>, nk<sub>max</sub>]</td><td>do</td></tr><tr><td>2</td><td></td><td>clusters, centers = agglom clustering(data)</td><td></td></tr><tr><td>3</td><td></td><td>div_cl = split clusters(clusters)</td><td></td></tr><tr><td>4</td><td></td><td>final_cl, finalCenters = merge clusters(div_cl)</td><td></td></tr><tr><td>5</td><td></td><td>end</td><td></td></tr></table>

## A. Reorganization using initial and final points

Two different approaches were taken to evaluate the best solution to further divide the clusters according to their initial and final points:

* Cluster both initial and final points on the same array, establishing automatically the sub-cluster groupings

* Cluster the initial point, then the final point and list the grouping created comparing both results.

Given that a search based on a number of clusters is already being executed, the mean-shift algorithm was used to execute these two options of post-processing, avoiding a nested search. After the merge process will happen, to fuse clusters with fairly similar characteristics, only differing from a few meters from each other while retaining its significant part, for example a left turn, a street crossing, etc.

Evaluating if two sub-clusters should be merged is done using the medoid of each cluster, together with the spread of elements in each sub-cluster, Equation (8). The variable $m_i$ represents the medoid of the cluster $C_i$ and $N_i$ its number of trajectories. To merge one sub-cluster with another it is necessary that the distance between both medoids be smaller than the sum of spreads for the respective clusters (line 5 of algorithm 2). To discard small differences, the medoid of one cluster is projected into the other (line 3 of algorithm 2), so that the similarity disregard any errors in tracking or even small differences in the start or terminus of the trajectory. If this condition is true, there is another to be fulfilled: the calculated projection should be equal or higher than a certain percentage of the original trajectory (which is defined in Table II for all cases examined here).

$$\text{spr}_{C_i} = \frac{1}{N_i} \sum_{x_i \in C_i} d_{\text{DTW}}(m_i, x_i) \quad (8)$$

Equations (9) and (10) show how one trajectory is projected onto another. For two generic trajectories $\mathbf{tr}_a = (p_{a0}, p_{a1}, \dots, p_{an})$ and $\mathbf{tr}_b = (p_{b0}, p_{b1}, \dots, p_{bn})$, two loops are executed: one for the initial point of $\mathbf{tr}_a$ and another its final point. Inside these loops the index $j$ in increased from 0 (or decreased from $N_b$ for the terminus) to find the interval of points where $p_{a0}$ (or $p_{an}$) is perpendicularly projected. The cutting point is obtained when $\lambda_{bj} \in [0, 1]$, meaning that the projection of $p_{a0}$ (or $p_{an}$), is between points $j$ and $j + 1$ from $\mathbf{tr}_b$. If no cutting point is detected then all the trajectory is used in the comparison.
<table><caption>**Algorithm 2:** Merging sub-clusters</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>for Cluster <i>i</i> in <i>C</i> do</td></tr><tr><td>2</td><td>for Cluster <i>j</i> ≠ <i>i</i> in <i>C</i> do</td></tr><tr><td>3</td><td>    proj<sub><i>i</i>→<i>j</i></sub> = Projection(medoid<sub><i>i</i></sub>, medoid<sub><i>j</i></sub>);</td></tr><tr><td>4</td><td>    d<sup><i>i</i>,<i>j</i></sup><sub>proj</sub> = DTW(med<sub><i>j</i></sub>, proj<sub><i>i</i>→<i>j</i></sub>);</td></tr><tr><td>5</td><td>    if d<sup><i>i</i>,<i>j</i></sup><sub>proj</sub> ≤ spr<sub><i>i</i></sub> + spr<sub><i>j</i></sub> then</td></tr><tr><td>6</td><td>        if trace<sub><i>i</i></sub> ≤ c<sub>min<sub><i>ac</i></sub></sub> · trace<sub><i>j</i></sub> then C<sub><i>j</i></sub> ← C<sub><i>i</i></sub>;</td></tr><tr><td>7</td><td>        else</td></tr><tr><td>8</td><td>    continue</td></tr></tbody></table>

$$v_{p_{a_0} \to p_{b_j}} = p_{b_j} - p_{a_0} \qquad (9)$$

$$\lambda_{b_j} = \frac{[v_{p_{a_0} \to p_{b_j}}] \cdot \hat{v}_{p_{b_j, j+1}}}{\|\hat{v}_{p_{b_j, j+1}}\|} \qquad (10)$$

## B. Evaluating cluster partition quality

Finally, it is necessary to evaluate the clusters according to the similarity of elements inside each cluster and in other clusters as well. Three metrics will be used for it: the Davies-Bouldin index the Silhouette score and the spread on cluster, proposed here. The DB index will be slightly modified to allow a better representation of the distribution quality for a certain number of clusters, while the latter will be used as it was defined in [17].

1) **Davies-Bouldin index (DB):** The Davies-Bouldin index (DB) is originally defined as the average of the maximal value of $R_{ij}$, as if defined in Equations (11) and (12). Equation (8) is used to calculate the spread $s_i$. Since it is the maximal value of $R_{ij}$ that is used to calculate the final score, it has creates a dependency to the number of clusters, i.e. the decrease on the score is connected to a higher number of cluster and not necessarily a better distribution.

$$R_{ij} = \frac{(s_i + s_j)}{d_{ij}} \qquad (11)$$

$$\text{DB}_{n_c} = \frac{1}{n_c} \sum_{i=0}^{n_c} \left[ \max_{j=1,...,n_c, i \neq j} R_{ij} \right] \qquad (12)$$

Thus, a small modification was done, to use the average of $R_{ij}$, not its maximal value, as it can be seen in Equation (13). This offers less bias to the number of clusters, considering always all the distribution of elements being evaluated. The $n_c - 1$ discounts the distance of the medoid to itself, which is zero.

$$\text{DB}_{n_c} = \frac{1}{n_c} \frac{1}{n_c - 1} \cdot \sum_{i=0}^{n_c} \sum_{j=0}^{n_c} R_{ij} \qquad (13)$$

2) **Silhouette score (Silh.):** Another metric to evaluate the clustering quality is the silhouette score, proposed in [17]. Differently from the DB score, it is calculated for each element being clustered, with $a(i)$ begin the average dissimilarity (DTW distance in this case) of element $i$ to all other

elements in its cluster. The other value necessary to calculate the silhouette is the minimum average dissimilarity between the element in question to the other clusters, Equation (15).

$$s_{C_i}(i) = \frac{b(i) - a(i)}{\max(a(i), b(i))} \qquad (14)$$

$$b(i) = \min_{C_j \neq C_i} d_{\text{DTW}}(b(i), C_j) \qquad (15)$$

With the $s(i)$ for each trajectory, the silhouette score for the clustering is the average score for all elements. This score is contained in $[-1, 1]$, with a score close to 1 being excellent (distances inside cluster are much smaller than distance between clusters). As the DB score it compares a infra-cluster spread measure with inter-cluster distances, but in a individual level. As it will be seen, in some situations where many different clusters co-exist close to each other (subsection IV-B), it will not be a representative measure of cluster quality.

3) **Spread on cluster (Spr.):** This measure is somewhat similar to Equation (9) but was changed to capture the biggest difference between two members of the same cluster. In Equation (16) the average of all the spreads divided by the number of members in the respective cluster define the metric.

$$\theta_C = \frac{1}{\|C\|} \cdot \sum_{i=0}^{\|C\|} \frac{\max_{j,k \in C_i} [d_{\text{DTW}}(x_j, x_k)]}{\|C_i\|} \qquad (16)$$

This metric will be specially important for the pedestrian case, where the silhouette score is not representative given the close proximity of multiple clusters.

# IV. RESULTS

## A. Methodology

Two different sources of data will be used to test the algorithms proposed here: the inD dataset [7] and the roundD dataset [8]. Both are obtained using a unnamed aerial vehicle (UAV), the former containing four different intersections and the latter three roundabouts, all in Germany. The information about the number of trajectories per scenario and per road user can be found in table I; the number of files refers to the data-batch file division for each scenario: in the inD case only the first three intersections were used for the VRU (for the vehicle case all intersection were tested) and for the roundD there are two files with two different roundabouts that are not used (only a few observations); the other scenarios are obtained from the observation of a third one (files 02 to 23). This data was split into three scenarios, given the number of trajectories (the clustering results of the three scenarios could be merged using the algorithm proposed in [2]). For the clustering execution all three road user trajectories from inD were used (in blue at Table I), while for the roundD only the vehicles' trajectories could be used, due to the low number of observations for pedestrians and cyclists (in red at Table I).

Table II shows the parameters used for the tests that will be presented next. The number of cluster established the interval
<table><caption>TABLE I: Number of trajectories per scenario per road user</caption><thead><tr><th>Scenario</th><th>Cars</th><th>Pedestrians</th><th>Cyclists</th></tr></thead><tbody><tr><td colspan="4" style="text-align: center;">inD Dataset</td></tr><tr><td>Sce. 0 - Files 00 to 06</td><td>1826</td><td>144</td><td>83</td></tr><tr><td>Sce. 1 - Files 07 to 17</td><td>2337</td><td>770</td><td>420</td></tr><tr><td>Sce. 2 - Files 18 to 29</td><td>2133</td><td>2005</td><td>1669</td></tr><tr><td>Sce. 3 - Files 30 to 32</td><td>1098</td><td>38</td><td>39</td></tr><tr><td colspan="4" style="text-align: center;">rounD Dataset</td></tr><tr><td>Sce. 0 - Files 02 to 08</td><td>3819</td><td>4</td><td>20</td></tr><tr><td>Sce. 1 - Files 09 to 15</td><td>2984</td><td>5</td><td>21</td></tr><tr><td>Sce. 2 - Files 16 to 23</td><td>3983</td><td>2</td><td>32</td></tr></tbody></table>

clusters tested by the methods, while the *Min. trace* refers to the minimal percentage that the projected medoid must have to be merged with another cluster (subsection III-A, algorithm 2).

<table><caption>TABLE II: Parameters used in clustering</caption><thead><tr><th>User</th><th>Scenario</th><th>Nb. clusters</th><th>Min. trace</th></tr></thead><tbody><tr><td rowspan="3">Pedestrian</td><td>inD sce. 0</td><td>[5, 32]</td><td>0.6</td></tr><tr><td>inD sce. 1</td><td>[15, 45]</td><td>0.6</td></tr><tr><td>inD sce. 2</td><td>[70, 110]</td><td>0.7</td></tr><tr><td rowspan="3">Cyclist</td><td>inD sce. 0</td><td>[5, 20]</td><td>0.6</td></tr><tr><td>inD sce. 1</td><td>[15, 45]</td><td>0.6</td></tr><tr><td>inD sce. 2</td><td>[20, 60]</td><td>0.6</td></tr><tr><td rowspan="5">Car</td><td>inD sce. 0</td><td>[5, 15]</td><td>0.67</td></tr><tr><td>inD sce. 1</td><td>[5, 20]</td><td>0.6</td></tr><tr><td>inD sce. 2</td><td>[5, 28]</td><td>0.4</td></tr><tr><td>inD sce. 3</td><td>[5, 22]</td><td>0.45</td></tr><tr><td>rounD sce. 0-2</td><td>[10, 25]</td><td>0.6</td></tr></tbody></table>

The choice for this specific datasets was motivated because both of them are captured by drone, not by an automated vehicle in the environment, which could modify the behaviors observed and also because the metadata present in the dataset allows the trajectories to be plotted in a realistic background image. But the method presented here could be used in any other ensemble of trajectories. Also, all the data is used as-is, no trajectory in the dataset is discarded beforehand only a normalization is done on each dimension of the trajectory before calculating the dissimilarity matrix (Equation 1). All the methods were implemented in Python.

## B. *Pedestrians*

There are three intersection scenarios for pedestrians, all in the inD dataset. The most important problem with pedestrians is the high variability of maneuvers, because it has an almost constraint-free environment to evolve and also due to detection and tracking errors during the dataset acquisition and post-processing. It is a real challenge for the clustering process to treat all these problems and produce a compact cluster set. The results for scenarios 0 and 2 can be seen in table III and table IV.

The abbreviation *Agglo* refers to the pure agglomerative clustering, A2MS to the agglomerative followed by two mean-shifts, on the initial and final points separately and A1MS to one mean-shift on the initial and final points on the same array. For all tables, III through XI, the time indicated

is the average clustering time per cluster calculated. In the column best $n_k$ the first value is the nominal value of the number of clusters used, and in parenthesis it is the final number of clusters, and in the last column the percentage of the original trajectories present in the final classification (in parenthesis the number of rejected trajectories).

Both scenarios are two extremes for the clustering process: one has few trajectories and not many options for a pedestrian to evolve in the environment and the other has much more of both. It is exactly of possible destinations and the possibility to use the same space in both directions that make the silhouette measure to fail when choosing the best method. Since the number of optimal clusters for each method is different, the spread on cluster, defined in III-B.3 is the best choice of criteria to select the best method overall and in the current case it indicates that the agglomerative clustering with two mean-shift applications (A2MS) is the best option for the scenario 0 and 2 (for the 1 as well, for space limitations it will not be shown here).

<table><caption>TABLE III: Results for scenario 0 of inD dataset</caption><thead><tr><th>Method</th><th>Avg. time</th><th>DB</th><th>Slh.</th><th>Spr.</th><th>Best n<sub>k</sub></th><th>Nb trajs.</th></tr></thead><tbody><tr><td>Agglo.</td><td>0.001 s</td><td>0.3164</td><td><strong>0.7012</strong></td><td>1.7561</td><td>32 (22)</td><td>93.06% (10)</td></tr><tr><td>A2MS</td><td>13.38 s</td><td><strong>0.2546</strong></td><td>0.5753</td><td><strong>1.4250</strong></td><td>23 (19)</td><td>79.86% (29)</td></tr><tr><td>A1MS</td><td>10.36 s</td><td>0.2837</td><td>0.5823</td><td>1.5823</td><td>19 (18)</td><td>79.86% (29)</td></tr><tr><td>PAM</td><td>16.07 s</td><td>0.5407</td><td>0.5810</td><td>3.2508</td><td>26 (24)</td><td>98.61% (2)</td></tr><tr><td>Dissi.</td><td>16.68 s</td><td>0.4820</td><td>0.4921</td><td>3.5506</td><td>32 (29)</td><td>97.92% (3)</td></tr></tbody></table>

For both scenarios the same observation can be made: the A2MS method has the lowest DB score and spread while the pure agglomerative method has the biggest silhouette score. This is exactly because the latter method some clusters that should be separated end up together while for the former one they usually can be separated by the second clustering. In both cases, clusters with a single element are discarded from the final distribution (and not accounted for during the calculation of the scores presented in each table). The exactly same observations can be made for the scenario 2².

<table><caption>TABLE IV: Results for scenario 2 of inD dataset</caption><thead><tr><th>Method</th><th>Time</th><th>DB</th><th>Slh.</th><th>Spr.</th><th>Best n<sub>k</sub></th><th>Nb trajs.</th></tr></thead><tbody><tr><td>Agglo.</td><td>0.04 s</td><td>0.4751</td><td><strong>0.5330</strong></td><td>1.7863</td><td>108 (77)</td><td>98.45% (31)</td></tr><tr><td>A2MS</td><td>57.75 s</td><td><strong>0.3500</strong></td><td>0.4693</td><td><strong>1.3265</strong></td><td>81 (77)</td><td>94.96% (101)</td></tr><tr><td>A1MS</td><td>34.90 s</td><td>0.4257</td><td>0.4610</td><td>1.6738</td><td>72 (77)</td><td>96.71% (66)</td></tr><tr><td>PAM</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Dissi.</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

Figure 3 shows an example of the effect of clustering with initial and final points to then merge the most similar trajectories later can have. In Figure 3a the result of the pure agglomerative cluster separated cluster 2 can be seen, where it mix an outlier (the white cross in the middle of the trajectory). The A2MS method (A1MS resulted in the same result) not only was able to separate the outliers (figure 3b) but also merged another pertinent cluster with it.

<footnote n="2">²For lack of time to process the data, the results for the PAM and dissi. method are not available at submission time; it does not change the final conclusion. At the final submission the values will be updated.</footnote>
<figure><img src="image_7.png" /><figcaption>(a) Agglo. method clust. 2</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) A2MS method cluster 0</figcaption></figure>

<caption>Fig. 3: Differences between methods for scenario 0</caption>

Differences between pure agglomerative against both mean-shift post-processing one are clearly visible, but in scenario 2 the results are more similar. This is probably due to the high number of samples, with helped the pure agglomerative method to sift through the outliers, but as it can be seen in Figure 4, not enough; it can also be seen the efficacy of the A2MS method to remove mixed clusters.

<figure><img src="image_9.png" /><figcaption>(a) Agglo. cluster 20 (b) A2MS. cluster 13</figcaption></figure>

<caption>Fig. 4: Differences between methods for scenario 2</caption>

Concerning the comparison with the PAM and dissimilarity methods, one can see that they are inefficient in both fronts being evaluated here: take more time to calculate and do not produce tight clusters, specially because of the outliers present in the scene. In [2], these methods were used to discover the different maneuvers of vehicles, but it must be highlighted that car's behaviors are much more constrained than pedestrians (and thus prone to outliers) and that the few outliers observed were removed before execution.

### C. Cyclists

The cyclist data was acquired at the same intersections than the pedestrians. From the set of trajectories given by the dataset the approximate cyclist behavior can be considered as somewhat between cars and pedestrians, with a constrained movement, but still able to access multiple parts of the road environment. For scenario 1 results (Table V) the main distinction that can be made from the pure agglo. and its modification is the opposite of what was observed with pedestrians.

<table><caption>TABLE V: Results for inD dataset scenario 1</caption><thead><tr><th>Method</th><th>Time</th><th>DB</th><th>Slh.</th><th>Spr.</th><th>Best n<sub>k</sub></th><th>Nb trajs.</th></tr></thead><tbody><tr><td>Agglo.</td><td>0.0036 s</td><td>0.6452</td><td>0.7746</td><td>0.5748</td><td>28 (12)</td><td>96.19% (16)</td></tr><tr><td>A2MS</td><td>5.44 s</td><td><strong>0.4818</strong></td><td><strong>0.8139</strong></td><td><strong>0.4854</strong></td><td>26 (12)</td><td>94.52% (23)</td></tr><tr><td>A1MS</td><td>3.86 s</td><td>0.5199</td><td>0.8090</td><td>0.5365</td><td>26 (12)</td><td>95.95% (17)</td></tr><tr><td>PAM</td><td>44.50 s</td><td>0.6040</td><td>0.0794</td><td>0.7684</td><td>38 (35)</td><td>99.29% (3)</td></tr><tr><td>Dissi.</td><td>60.80 s</td><td>0.6666</td><td>0.1560</td><td>0.8571</td><td>42 (34)</td><td>98.10% (8)</td></tr></tbody></table>

Some of the trajectories were split between multiple clusters for the pure agglomerative method while for the A\*MS (meaning both A1MS and A2MS) methods these clusters

could be merged together, specially in scenario 2. There are other instances of this same behavior in different clusters as well. Both A\*MS methods can attribute their superior scores to the ability to remove outliers from clusters, as illustrated in Figure 5.

<table><caption>TABLE VI: Results for inD dataset scenario 2</caption><thead><tr><th>Method</th><th>Time</th><th>DB</th><th>Slh.</th><th>Spr.</th><th>Best n<sub>k</sub></th><th>Nb trajs.</th></tr></thead><tbody><tr><td>Agglo.</td><td>0.035 s</td><td>0.7789</td><td><strong>0.7286</strong></td><td>1.6400</td><td>54 (40)</td><td>99.16% (14)</td></tr><tr><td>A2MS</td><td>18.50 s</td><td><strong>0.5761</strong></td><td>0.6582</td><td><strong>1.0211</strong></td><td>70 (37)</td><td>95.21% (80)</td></tr><tr><td>A1MS</td><td>11.46 s</td><td>0.6484</td><td>0.6316</td><td>1.3435</td><td>61 (48)</td><td>97.84% (36)</td></tr><tr><td>PAM</td><td>1121.17 s</td><td>0.7331</td><td>0.2814</td><td>1.2032</td><td>60 (59)</td><td>99.94% (1)</td></tr><tr><td>Dissi.</td><td>983.17 s</td><td>0.8621</td><td>0.3115</td><td>1.3177</td><td>59 (56)</td><td>99.82% (3)</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>(a) Agglo. cluster 0 (b) A2MS cluster 19</figcaption></figure>

<caption>Fig. 5: Clustering differences for scenario 2</caption>

### D. Vehicles

For vehicles the volume of data increases, with the addition of inD dataset scenario 3 and the entire rounD dataset. Since the movements for vehicles are very constrained there are almost none eccentric behavior, hence the goal here is to eliminate all the erroneous samples; for example, trajectories that end at the middle of the intersection. In some cases this was possible, notably on scenarios 1 and 2 for the inD dataset, however, in scenario 0 one maneuver got separated as the result of the mean-shift and merge mechanism for the A2MS method (Figure 6). Beyond that all other maneuvers for 0 were correctly determined.

<table><caption>TABLE VII: Results for inD dataset scenario 0</caption><thead><tr><th>Method</th><th>Time</th><th>DB</th><th>Slh.</th><th>Spr.</th><th>Best n<sub>k</sub></th><th>Nb trajs.</th></tr></thead><tbody><tr><td>Agglo.</td><td>0.042 s</td><td>1.6403</td><td><strong>0.9014</strong></td><td>0.4538</td><td>15 (9)</td><td>99.67% (6)</td></tr><tr><td>A2MS</td><td>6.22 s</td><td><strong>0.9366</strong></td><td>0.8615</td><td><strong>0.1850</strong></td><td>23 (10)</td><td>99.23% (14)</td></tr><tr><td>A1MS</td><td>3.67 s</td><td>1.1694</td><td>0.6796</td><td>0.3095</td><td>22 (13)</td><td>99.45% (10)</td></tr><tr><td>PAM</td><td>50.07 s</td><td>2.5422</td><td>0.6283</td><td>0.6289</td><td>13 (13)</td><td>100% (0)</td></tr><tr><td>Dissi.</td><td>51.69 s</td><td>1.9416</td><td>0.3239</td><td>0.4938</td><td>15 (12)</td><td>99.84 (3)%</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>(a) A2MS method cluster 3 (b) A2MS method cluster 5</figcaption></figure>

Fig. 6: Problems with the A2MS merge when removing outliers

For scenario 1 the agglo., A2MS and A1MS were spot on, with the sole difference that two clusters detected by the agglo and the A1MS are actually outliers and were rejected by A2MS before they formed clusters.
<table><caption>TABLE VIII: Results for inD dataset scenario 1</caption><thead><tr><th>Method</th><th>Time</th><th>DB</th><th>Slh.</th><th>Spr.</th><th>Best n<sub>k</sub></th><th>Nb trajs.</th></tr></thead><tbody><tr><td>Agglo.</td><td>0.10 s</td><td>1.3758</td><td>0.8942</td><td>0.5533</td><td>19 (14)</td><td>99.79% (5)</td></tr><tr><td>A2MS</td><td>7.74 s</td><td><strong>0.8631</strong></td><td><strong>0.9044</strong></td><td><strong>0.2221</strong></td><td>20 (12)</td><td>99.23% (28)</td></tr><tr><td>A1MS</td><td>4.60 s</td><td>1.2074</td><td>0.8260</td><td>0.4431</td><td>18 (14)</td><td>99.19% (29)</td></tr><tr><td>PAM</td><td>134.08 s</td><td>1.5372</td><td>0.7375</td><td>0.4658</td><td>19</td><td>100.0%</td></tr><tr><td>Dissi.</td><td>149.07 s</td><td>2.2906</td><td>0.7046</td><td>0.5292</td><td>20</td><td>100.0%</td></tr></tbody></table>

As for the PAM and dissi methods, they split different maneuvers into different clusters: Figure 7a represents the cluster that was divided into Figures 7c, 7d, 7e, which explains huge the disparity shown in Table VIII.

<figure><img src="image_12.png" /><figcaption>(a) A2MS method cluster 0</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) A2MS method cluster 3</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c) PAM cluster 12</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d) PAM cluster 16</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(e) PAM cluster 9</figcaption></figure>

<caption>Fig. 7: Outliers (Agglo.) and clusters from A2MS and PAM</caption>

Table IX marks is the first time that the A1MS method had a better DB index and its analogue, due to the splitting of a curve maneuver that contained many samples. This is why the DB score is not used to define the best method, even when the number of cores is the same: there is some situations to split a cluster that should be an unit might be beneficial because of the spread (Equation (11)) calculation. Besides that, again the pure agglomerative cluster is not capable to split maneuvers that share most of their length (Figures 8a and 8b)

<table><caption>TABLE IX: Results for inD dataset scenario 2</caption><thead><tr><th>Method</th><th>Time</th><th>DB</th><th>Slh.</th><th>Spr.</th><th>Best n<sub>k</sub></th><th>Nb trajs.</th></tr></thead><tbody><tr><td>Agglo.</td><td>0.056 s</td><td>1.7995</td><td>0.8425</td><td>1.3679</td><td>17 (13)</td><td>99.81% (4)</td></tr><tr><td>A2MS</td><td>8.02 s</td><td>1.1375</td><td><strong>0.8466</strong></td><td><strong>0.7153</strong></td><td>22 (19)</td><td>98.97% (22)</td></tr><tr><td>A1MS</td><td>5.06 s</td><td><strong>1.1286</strong></td><td>0.8135</td><td>0.8857</td><td>15 (19)</td><td>99.02% (21)</td></tr><tr><td>PAM</td><td>251.27 s</td><td>1.7625</td><td>0.3552</td><td>1.2958</td><td>24</td><td>100% (0)</td></tr><tr><td>Dissi.</td><td>209.87 s</td><td>1.9005</td><td>0.3444</td><td>1.5167</td><td>22 (21)</td><td>99.95% (1)</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>(a) Agglo. cluster 0</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(b) Agglo. cluster 2</figcaption></figure>

<caption>Fig. 8: Differences for scenario 2</caption>

For scenario 3 the A1MS actually is the better option, but

there is only a difference of two samples classified differently from this case.

<table><caption>TABLE X: Results for inD dataset scenario 3</caption><thead><tr><th>Method</th><th>Time</th><th>DB</th><th>Slh.</th><th>Spr.</th><th>Best n<sub>k</sub></th><th>Nb trajs.</th></tr></thead><tbody><tr><td>Agglo.</td><td>0.014 s</td><td><strong>1.9011</strong></td><td><strong>0.9285</strong></td><td>0.2949</td><td>18 (10)</td><td>99.27% (2)</td></tr><tr><td>A2MS</td><td>4.46 s</td><td>2.3577</td><td>0.8696</td><td>0.2314</td><td>20 (9)</td><td>99.18% (9)</td></tr><tr><td>A1MS</td><td>1.98 s</td><td>2.1899</td><td>0.8671</td><td><strong>0.2289</strong></td><td>19 (9)</td><td>99.18% (9)</td></tr><tr><td>PAM</td><td>237.22 s</td><td>3.1968</td><td>0.6639</td><td>0.5730</td><td>19 (12)</td><td>99.82% (2)</td></tr><tr><td>Dissi.</td><td>64.44 s</td><td>3.1112</td><td>0.6513</td><td>0.5730</td><td>18 (12)</td><td>99.82% (2)</td></tr></tbody></table>

Since the trajectories for the rounD dataset are fairly different in length and direction the DTW distance measure is able to really account trajectories from different maneuvers, as it can be seen in table XI.

<table><caption>TABLE XI: Results for rounD dataset scenario 0</caption><thead><tr><th>Method</th><th>Time</th><th>DB</th><th>Slh.</th><th>Spr.</th><th>Best n<sub>k</sub></th><th>Nb trajs.</th></tr></thead><tbody><tr><td>Agglo.</td><td>0.57 s</td><td>0.4237</td><td><strong>0.9470</strong></td><td>0.1185</td><td>22 (14)</td><td>99.79% (8)</td></tr><tr><td>Agglo. 2MS</td><td>10.45 s</td><td><strong>0.2817</strong></td><td>0.6721</td><td><strong>0.1135</strong></td><td>19 (24)</td><td>99.71% (11)</td></tr><tr><td>Agglo. MS</td><td>5.72 s</td><td>0.3223</td><td>0.7633</td><td>0.1224</td><td>22 (22)</td><td>99.63% (14)</td></tr><tr><td>PAM</td><td>329.96 s</td><td>0.4174</td><td>0.5618</td><td>0.5193</td><td>25 (24)</td><td>99.97% (1)</td></tr><tr><td>Dissi.</td><td>424.83 s</td><td>0.6701</td><td>0.7870</td><td>0.1184</td><td>24 (16)</td><td>99.97% (1)</td></tr></tbody></table>

But there is something that are not as salient in the short turns in the intersection from inD dataset but is in this case. Given the size of the roundabout the position in which the vehicles execute the trajectory becomes a discriminating parameter, i.e. the clusters also accounted if they are on the inside or outside (Figures 9a and 9b comparing with 9c) together with the lane in which the vehicle ends or starts its trajectory. This difference mostly impacted the merge step given that now there is a lateral distance through the curve that is bigger than the spread over both clusters. If this division is appropriate or not is in the eye of the beholder.

<figure><img src="image_19.png" /><figcaption>(a) Agglo. cluster 0 (b) A2MS cluster 19 (c) A2MS cluster 32</figcaption></figure>

<caption>Fig. 9: Clustering differences for scenario 1</caption>

As a general comment, for the current use-case, the method A2MS proved vastly better than any other, which was made clear by the spread on cluster score defined here. It captures the tightness of each cluster much better than the DB index, which ends up being translated as clusters with few to none outliers in their midst.

V. CONCLUSION

A new method to cluster trajectories, A2MS, together with a metric defined for the trajectory clustering case, the spread on cluster, were proposed and tested with the datasets inD and rounD. Using an hierarchical clustering combined with DTW distance measure and a as cluster distribution measure, A2MS proved to be the most efficient in the tested methods
to produce tight and concentrated clusters with minimal number of outliers.

The immediate next step is to use this clustering method in conjunction with the longitudinal approach proposed by [2] to extract drivers' behaviors from real data. But more broadly the method proposed here has multiple uses, from prepare data to learning tasks for planning, decision-making or prediction to even the study of traffic flow in a predetermined zone. Ultimately, this method allows in the future to collect data to train a representation of trajectories so that comparisons could be made with trajectories from different road configurations.

## REFERENCES

[1] M. Y. Choong, L. Angeline, R. K. Yin Chin, K. Beng Yeo, and K. T. Kin Teo, "Modeling of vehicle trajectory clustering based on less for traffic pattern extraction," in *2017 IEEE 2nd International Conference on Automatic Control and Intelligent Systems (I2CACIS)*. IEEE, 2017, pp. 74-79.

[2] N. de Moura and F. Nashashibi, "Extraction of vehicle behaviors at intersections," in *IEEE 2023 Intelligent Transportation Systems Conference*, Bilbao, Spain, Sep. 2023. [Online]. Available: https://inria.hal.science/hal-04314028

[3] W. Wang, F. Xia, H. Nie, Z. Chen, Z. Gong, X. Kong, and W. Wei, "Vehicle trajectory clustering based on dynamic representation learning of internet of vehicles," *IEEE Transactions on Intelligent Transportation Systems*, vol. 22, no. 6, pp. 3567-3576, 2021.

[4] L. Ries, P. Rigoll, T. Braun, T. Schulik, J. Daube, and E. Sax, "Trajectory-based clustering of real-world urban driving sequences with multiple traffic objects," in *2021 IEEE International Intelligent Transportation Systems Conference (ITSC)*, 2021, pp. 1251-1258.

[5] M. Y. Choong, L. Angeline, R. K. Y. Chin, K. B. Yeo, and K. T. K. Teo, "Modeling of vehicle trajectory using k-means and fuzzy c-means clustering," in *2018 IEEE International Conference on Artificial Intelligence in Engineering and Technology (IICAIET)*. IEEE, 2018, pp. 1-6.

[6] H. Liu, J. Li, Y. Wu, and Y. Fu, "Clustering with outlier removal," *IEEE Transactions on Knowledge and Data Engineering*, vol. 33, no. 6, pp. 2369-2379, 2021.

[7] J. Bock, R. Krajewski, T. Moers, S. Runde, L. Vater, and L. Eckstein, "The ind dataset: A drone dataset of naturalistic road user trajectories at german intersections," in *2020 IEEE Intelligent Vehicles Symposium (IV)*. IEEE, 2020, pp. 1929-1934.

[8] R. Krajewski, T. Moers, J. Bock, L. Vater, and L. Eckstein, "The round dataset: A drone dataset of road user trajectories at roundabouts in Germany," in *2020 IEEE 23rd International Conference on Intelligent Transportation Systems (ITSC)*, 2020, pp. 1-6.

[9] E. Barmpounakis and N. Geroliminis, "On the new era of urban traffic monitoring with massive drone data: The pneuma large-scale field experiment," *Transportation Research Part C: Emerging Technologies*, vol. 111, pp. 50-71, 2020. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0968090X19310320

[10] W. Zhan, L. Sun, D. Wang, H. Shi, A. Clausse, M. Naumann, J. Kummerle, H. Konigshof, C. Stiller, A. de La Fortelle et al., "Interaction dataset: An international, adversarial and cooperative motion dataset in interactive driving scenarios with semantic maps," *arXiv preprint arXiv:1910.03088*, 2019.

[11] M. Müller, *Information retrieval for music and motion*. Springer, 2007, vol. 2.

[12] H. Sakoe and S. Chiba, "Dynamic programming algorithm optimization for spoken word recognition," *IEEE Transactions on Acoustics, Speech, and Signal Processing*, vol. 26, no. 1, pp. 43-49, 1978.

[13] F. Itakura, "Minimum prediction residual principle applied to speech recognition," *IEEE Transactions on Acoustics, Speech, and Signal Processing*, vol. 23, no. 1, pp. 67-72, 1975.

[14] F. Petitjean, A. Ketterlin, and P. Gançarski, "A global averaging method for dynamic time warping, with applications to clustering," *Pattern Recognition*, vol. 44, no. 3, pp. 678-693, 2011. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S003132031000453X

[15] H. Schütze, C. D. Manning, and P. Raghavan, *Introduction to information retrieval*. Cambridge University Press Cambridge, 2008, vol. 39.

[16] L. Kaufman and P. J. Rousseeuw, *Finding groups in data: an introduction to cluster analysis*. John Wiley & Sons, 2009.

[17] P. J. Rousseeuw, "Silhouettes: A graphical aid to the interpretation and validation of cluster analysis," *Journal of Computational and Applied Mathematics*, vol. 20, pp. 53-65, 1987. [Online]. Available: https://www.sciencedirect.com/science/article/pii/0377042787901257