Maël Primet
email: mael.primet@parisdescartes.fr

Lionel Moisan
email: lionel.moisan@parisdescartes.fr

Point tracking: an a-contrario approach

Keywords: point tracking, a-contrario detection, motion correspondence 1

In this work, we propose a new approach to recover trajectories from points previously detected in a sequence of images. In presence of spurious and missing detections, actual trajectories can be characterized by an acontrario model, that estimates the probability of observing a similar trajectory in random data. This results in a single criterion combining trajectory characteristics (duration, number of points, smoothness) and data statistics (number of images and detected points), which can then be used to drive a dynamic programming algorithm able to extract sequentially the most meaningful trajectories. The performances obtained on synthetic and real-world data are studied in detail, and shown to compare favorably to the state-of-the-art ROADS algorithm.

in accelerators [START_REF] Cornelissen | The new atlas track reconstruction (newt)[END_REF]. Numerous variations on this problem have been formulated, and several algorithms have been developed to try to solve them. A common strategy (see [START_REF] Yilmaz | Object tracking: A survey[END_REF] for a recent review) is to detect the objects in each image and associate them with a geometric representation, such as points (Veenman et al, 2003b;[START_REF] Serby | Probabilistic object tracking using multiple features[END_REF], geometric shapes [START_REF] Comaniciu | Kernel-based object tracking[END_REF], outlines [START_REF] Yilmaz | Contour-Based Object Tracking with Occlusion Handling in Video Acquired Using Mobile Cameras[END_REF], skeletal models [START_REF] Ali | Segmentation and recognition of continuous human activity[END_REF], and with appearance features decribed, e.g., by templates [START_REF] Fieguth | Color-Based Tracking of Heads and Other Mobile Objects at Video Frame Rates[END_REF], active appearance models [START_REF] Edwards | Interpreting face images using active appearance models[END_REF], or probability densities of object appearance [START_REF] Zhu | Region Competition: Unifying Snakes, Region Growing, and Bayes/MDL for Multiband Image Segmentation[END_REF][START_REF] Paragios | Geodesic Active Contours and Level Sets for the Detection and Tracking of Moving Objects[END_REF]. The detected objects are then tracked across frames using an algorithm that tightly depends on the object representation. According to [START_REF] Yilmaz | Object tracking: A survey[END_REF], tracking algorithms can be broadly classified in three categories:

1. Point tracking (Veenman et al, 2003b;[START_REF] Reid | An algorithm for tracking multiple targets[END_REF][START_REF] Bar-Shalom | Sonar tracking of multiple targets using joint probabilistic data association[END_REF][START_REF] Streit | Maximum likelihood method for probabilistic multihypothesis tracking[END_REF][START_REF] Shafique | A non-iterative greedy algorithm for multi-frame point correspondence[END_REF]. Objects are represented as points and are generally tracked across frames by evolving their state (object position and motion); 2. Kernel tracking [START_REF] Comaniciu | Kernel-based object tracking[END_REF][START_REF] Tao | Object Tracking with Bayesian Estimation of Dynamic Layer Representations[END_REF][START_REF] Black | EigenTracking: Robust Matching and Tracking of Articulated Objects Using a View-Based Representation[END_REF][START_REF] Avidan | Support vector tracking[END_REF]. Objects are represented by a combination of shape and appearance, for instance an ellipse with a color histogram. They are then tracked by computing the motion of the kernel in consecutive frames, often modeled with parametric transforms such as translations and rotations; 3. Silhouette tracking [START_REF] Blake | Active Contours: The Application of Techniques from Graphics,Vision,Control Theory and Statistics to Visual Tracking of Shapes in Motion, 1st edn[END_REF][START_REF] Bertalmío | Morphing active contours[END_REF][START_REF] Ronfard | Region-based strategies for active contour models[END_REF][START_REF] Kang | Object reacquisition using invariant appearance model[END_REF][START_REF] Huttenlocher | Tracking nonrigid objects in complex scenes[END_REF][START_REF] Sato | Temporal spatio-velocity transform and its application to tracking and interaction[END_REF]. Objects regions are estimated in each frame, and are usually tracked by either shape matching or contour evolution.

In this work, we restrict ourselves to the tracking of objects as points, without appearance information. We assume that these points have already been detected in the sequence, but imperfectly, in the sense that we may have to deal with both spurious points and missing detections. Given such a sequence of frames containing the detected points in each image, the goal is to extract the trajectories as lists of points appearing in successive frames, possibly separated by holes (missing detections), while avoiding spurious points. As usual, we will assume that a given point can belong to only one trajectory, which means that point collisions are ruled by an occlusion principle.

Related work

The most well-known point tracking algorithm certainly is the Multiple Hypothesis Tracker (MHT) of [START_REF] Reid | An algorithm for tracking multiple targets[END_REF], which, in theory, enumerates all possible trajectory combinations of the observed points, and selects the one having the maximal likelihood (a probabilistic motion model being given). This problem is indeed NP-hard, and leads to exponentially many trajectory combination lookups (if there are n frames and m points, and we assume there is no occlusion, noise points, or objects leaving or entering the scene, there are already m! n-1 possible trajectory combinations), and thus approximate solutions and heuristics are needed to accelerate the search (often requiring thresholds to prune the search tree early). Moreover, this complex model implies parameter tuning to optimize the underlying motion model and the efficiency/coverage tradeoff of the heuristics.

To overcome the exponential growth of the state space, researchers have proposed a wealth of heuristics and approximations to the point tracking problem over the years. Such an approximation is the Joint Probabilistic Data Association Filter (JPDAF) proposed by [START_REF] Bar-Shalom | Sonar tracking of multiple targets using joint probabilistic data association[END_REF], which relaxes the hypothesis that the points must be disjoint. Instead of assigning each track ending in frame k -1 to a particular object in frame k as in MHT (and thus having several possible hypotheses, resulting in an increase of complexity), the JPDAF algorithm assigns to every track ending in frame k -1 a weighted combination of all points of frame k, depending on a likelihood estimate with respect to a predicted position. However, this approach assumes that the number of objects tracked in the images is constant, and the relaxation of the disjointness hypothesis leads to trajectory mergings, which is often an undesired feature. [START_REF] Sethi | Finding trajectories of feature points in a monocular image sequence[END_REF] propose to solve the correspondence problem greedily. They initialize the trajectories using the nearest-neighbor criterion, and then improve the current solution by exchanging correspondences between frames in order to minimize the global cost. They also propose a modified algorithm that alternates between forward and backward passes through the sequence to help mitigate the problem of the nearest neighbor initialization. Their approach is much faster than MHT, but does not permit to take noise, occlusions, entries or exits into account. [START_REF] Salari | Feature point correspondence in the presence of occlusion[END_REF] address some shortcomings of the previous method, namely the fact that it assumes a constant number of points in the sequence, and that there can be no entries or exits of objects in the scene. It therefore allows for occlusion, entry and exit of points, as well as the presence of spurious points. Each trajectory is made of either points detected in a frame, or "phantom points", which correspond to added (interpolated) points when there are missing detections. Note that their approach involves two parameters (the maximum allowed speed and the maximum allowed acceleration) which may be difficult to set. [START_REF] Rangarajan | Establishing motion correspondence[END_REF] propose to solve the problem by using a proximal uniformity constraint which combines requirements on the maximum speed and the acceleration of objects. They propose to make the assigmnent choices in an order driven by a notion of minimal risk, still in a greedy way. They use an optical flow algorithm to initialize the point motion between the first two images, and deal with occlusion and missing detections by using linear interpolation. They do not allow for spurious points, or objects leaving or entering the scene. Veenman, Reinders, and Backer (2003b) propose a greedy tracking algorithm called ROADS, which is capable of handling missing and spurious detections, as well as entries and exits of objects. Rather than optimizing a global cost, they consider a restricted temporal scope (usually two or three frames forward), and find the optimal assignments minimizing the cost on these frames. Since the restricted problem is still NP-hard, they have some heuristics that help them prune the search tree. They keep the assignment between the two first frames of the local scope, and iterate the process on the following frames. The assignment between the first two frames of the sequence is initialized using the nearest neighbor criterion, and the effect if this approximation is mitigated as in [START_REF] Sethi | Finding trajectories of feature points in a monocular image sequence[END_REF]) by a forward and a backward pass. The ROADS algorithm is an extension of the previous well-established GOA algorithm [START_REF] Veenman | Resolving motion correspondence for densely moving points[END_REF], which itself compares favorably to three classical algorithms [START_REF] Salari | Feature point correspondence in the presence of occlusion[END_REF], [START_REF] Rangarajan | Establishing motion correspondence[END_REF] and [START_REF] Chetverikov | Feature point tracking for incomplete trajectories[END_REF]). The ROADS algorithm has been shown to outperform both GOA and the MHT algorithm of [START_REF] Reid | An algorithm for tracking multiple targets[END_REF].

Another state-of-the-art algorithm, developped for fluorescence particle tracking, is described in [START_REF] Sbalzarini | Feature point tracking and trajectory analysis for video imaging in cell biology[END_REF]. However, it addresses a different (and easier) tracking problem, since it requires information on the intensity of the detections. Moreover, it uses a prior motion model (nearest neighbor) that is more adapted to small displacements than to the kind of smooth trajectories (low acceleration) that we will consider in the present paper. [START_REF] Fleuret | Multi-Camera People Tracking with a Probabilistic Occupancy Map[END_REF] propose to track multiple persons in multiple camera views using a probabilistic map of the individuals locations, coupled with a Dynamic Programming algorithm that tracks each person in isolation, rather than conjointly. They use both an appearance model and a motion model to describe the objects to track.

In essence, the approaches above try to restrict the search space either by constraining the optimization to local choices or by limiting the simultaneous number of objects tracked. A recent and promising approach introduced in [START_REF] Jiang | A linear programming approach for multiple object tracking[END_REF] tries to solve the original unconstrained problem by optimizing a global criterion using Linear Programming, where the correspondence decisions are not binary choices, but continuous values in [0, 1], rendering the optimization problem convex, and thus efficiently computable. In practice, the estimated values are almost always equal to either 0 or 1, and it is easy to convert them into disjoint trajectories. The algorithm assumes that the number of objects in the images is constant, but has been later extended in [START_REF] Berclaz | Multiple Object Tracking using Flow Linear Programming[END_REF]; [START_REF] Berclaz | Multiple object tracking using k-shortest paths optimization[END_REF] to accommodate a variable number of objects that can enter and leave the scene in prespecified locations. The authors also prove that when there exists a unique global optimum, it is necessarily a boolean optimum.

Trajectory estimation versus trajectory detection

These algorithms have some common limitations. First, when they take a varying number of objects into account, as well as spurious and missing detections, they face the difficulty of choosing an appropriate global cost for their optimization problem. Indeed, they need a penalization for spurious points, and most of the time they will simply introduce a fixed cost for them. This creates a subtle (and quite arbitrary) interplay between the cost of spurious points and the cost of detected trajectories, which is not easy to control and grasp. Second, these algorithms often have many parameters, which is fine in theory, but quickly becomes a hassle when one needs to set them for each practical use. Last, as they all try to solve an optimization problem, these algorithms suffer from a classical flaw: they always find something, since they try to find the best explanation of the data in terms of some structure, without trying first to prove that the structure is present. In particular, all these algorithms will, for some values of the parameters, find trajectories in random data made of pure random points (without motion coherence).

In the present work, we propose a new approach for trajectory detection, which can guarantee that no trajectory will be found in general in such random data. This work is based on the a-contrario framework [START_REF] Desolneux | From Gestalt Theory to Image Analysis: A Probabilistic Approach[END_REF], which permits to derive absolute detection thresholds. These thresholds are then used to drive an algorithm that is able to analyze trajectories globally in time while avoiding the three aforementionned limitations.

In Section 2, we first consider trajectories without holes, that is, the case where no data point is missing (but spurious points are expected). After recalling the basic principles of the a-contrario statistical framework, we derive an explicit criterion for trajectory detection and present an algorithm based on Dynamic Programming [START_REF] Bellman | The theory of dynamic programming[END_REF]. We also analyze some theoretical consequences of the acontrario thresholds, in particular the link between the number of points, the number of images and the maximum allowed acceleration. Then Section 3 extends the theory and the algorithm of Section 2 to the more general case of trajectories that contain holes. In Section 4, the state-of-theart ROADS algorithm (Veenman et al, 2003b) is considered and various experiments (following, for most of them, the methodology proposed in the original ROADS paper) are led to compare its performances with the a-contrario algorithm we propose. Aside from a very convenient reduction of the number of parameters (1 for the NFA algorithm, versus 4 in the ROADS experiments), the a-contrario algorithm significantly outperforms the ROADS algorithm in terms of precision, robustness, and sensitivity to parameters, both in the "no-hole" and in the "hole" versions. Experiments are also conducted on real data, namely a snow sequence (that we make publicly available online) for which the ground truth has been manually obtained. Again, the results are clearly in favor of the a-contrario algorithm. We finally conclude in Section 5, and comment on the strengths, limitations, and perspectives offered by the present work.

Trajectories without holes

In this part, we consider trajectories without holes, that is, we assume that there are no missing detections (but possibly spurious detections, and points leaving and entering the scene).

Principles of the a-contrario framework

The trajectory detection method that we propose relies on the a-contrario framework introduced by Desolneux, Moisan and Morel (see [START_REF] Desolneux | From Gestalt Theory to Image Analysis: A Probabilistic Approach[END_REF] for a recent presentation). The idea underlying its development (dubbed "Helmholtz Principle") is that the human visual system detects structures in an image as coincidences which could not appear by chance in a random setting. Conceived at first to detect structures issued from Gestal Theory [START_REF] Wertheimer | Untersuchungen zur lehre von der gestalt[END_REF][START_REF] Kanizsa | Grammatica del vedere[END_REF], this methodology has been applied to a large variety of image processing tasks, aiming at detecting structures like alignments [START_REF] Desolneux | Meaningful alignments[END_REF], edges [START_REF] Desolneux | Edge detection by Helmholtz principle[END_REF], stereo coherence [START_REF] Moisan | A probabilitic criterion to detect rigid point matches between two images and estimate the fundamental matrix[END_REF], spots [START_REF] Grosjean | A-contrario detectability of spots in textured backgrounds[END_REF], image changes [START_REF] Robin | An A-Contrario Approach for Subpixel Change Detection in Satellite Imagery[END_REF], etc. Note that it has also been used in the context of Motion Analysis to cluster motioncoherent regions in image sequences [START_REF] Veit | Space-time a contrario clustering for detecting coherent motion[END_REF], or to distinguish moving regions from the scene background [START_REF] Veit | An a contrario decision framework for region-based motion detection[END_REF].

We here recall the formalization of the a-contrario framework as it was presented in [START_REF] Grosjean | A-contrario detectability of spots in textured backgrounds[END_REF]. The a-contrario methodology is based on two ingredients: a naive model, and one or several measurements defined on the structures of interest. The naive model describes typical situations where no structure should be detected. For instance, when trying to discover alignments of points in an image, a naive model could consist of uniform and independent draws of the point locations, where no interesting structure would usually appear (see Fig. 1).

Fig. 1 Illustration of Hemholtz principle. Why can't we help seeing an alignment of dots on the left image? According to Helmholtz principle, we a priori assume that the dots should have been drawn uniformely and independently as in the right image, and we perceive a structure (here a group of aligned dots) because such an alignment is very unlikely to happen by chance. Alignments of three dots can be found in the right image, but they do not pop out, because they are likely to happen by chance considering the total number of points. The formalization of this principle is realized in a-contrario detection models.

To detect structures (e.g. alignments of points) in data using Helmholtz Principle, we need to define in what way an observation can be significant. If the measurement function is high when the structure is pregnant, we can relate the "amount of surprise" when observing the measurement x to the probability P(X x), where X is the random variable corresponding to the distribution of x in the naive model. We will usually have several measurements (x i ) i∈I (in the example above, one for each possible alignment), and in the classical a-contrario framework the amount of surprise will be measured by a number of false alarms. More formally, we have the following Definition 1 (Number of False Alarms) Let (X i ) 1 i N be a set of random variables. A family of functions F i (x) i is a NFA (number of false alarms) for the random variables

(X i ) i if ∀ε > 0, E #{i, F i (X i ) ε} ε (1)
(as usual, the notation "#S" stands for the cardinal of the set S).

A measurement x i such that F i (x i ) ε is said to be detected at level ε, or ε-meaningful. We say that a measurement is meaningful if it is 1-meaningful. This number of false alarms ensures that the average number of detections made in the naive model (that is, false detections) at level ε is less than ε.

The classical way to construct a NFA is given by the following proposition (see [START_REF] Grosjean | A-contrario detectability of spots in textured backgrounds[END_REF]).

Proposition 1 (NFA construction) Let (X i ) 1 i N be a set of random variables and (w i ) 1 i N a set of positive real numbers, then the function

NFA(i, x i ) = w i • P(X i x i ) (2)
is a NFA as soon as

N ∑ i=1 1 w i 1 and in particular if w i = N for all i. Remark 1 (NFA approximation) If (F i ) i is a NFA, then any family of functions (G i ) i verifying F i (t) G i (t) for all t is still a NFA. Hence, a function satisfying NFA(i, x i ) w i • P(X i x i )
will define a NFA as soon as ∑ N i=1 1/w i 1.

Trajectory detection

We are given a sequence of K images, each containing N points (to ease the notations we consider a constant number of points throughout the sequence but everything can be smoothly extended to the non-constant case as will be shown later), and whose support domain is taken to be the square [0, 1] × [0, 1] (again, the method can be adapted to arbitrary image sizes, as shown in Section 2.3). Following Helmholtz Principle, the naive model will here be a random uniform draw of N points in each of the K images (intuitively, we should not see trajectories appearing in the realizations of this model). The associated i.i.d. uniformely distributed random variables corresponding to the points of image k (1 k K) will be denoted by X k 1 , ..., X k N . We now define the structures of interest.

Definition 2 (Trajectories without holes) A trajectory T of length ℓ starting at frame k 0 is a tuple T = (k 0 , i 1 , ..., i ℓ ), where 1 i p N for all p and 1 ℓ Kk 0 + 1. We will denote by T the set of all trajectories. There is a natural equivalence between a trajectory T ∈ T and the tuple of variables X T = (X

k 0 i 1 , ..., X k 0 +ℓ-1 i ℓ
) that we shall therefore sometimes abusively call a (random) trajectory too.

A realization t of the random variable X T will be called a realization of the trajectory T . We have to keep in mind that we are working with the naive model (where points are randomly distributed), and thus, a realization of X T should not look like what we would intuitively call a trajectory.

The second ingredient of the a-contrario model is the measurement function. For most applications, a natural choice is to define a "good trajectory" as a smooth one, and more precisely as a trajectory with a small acceleration. The proposed model could as well be easily adapted to a variety of other measurements; for instance, we could define a "good trajectory" as one that moves slowly, and take the measurement function to be the maximal speed of the object along the trajectory. Let us stick to the idea of good trajectories being those bearing a small acceleration. We still have many ways to define this, the two most obvious choices being to control the maximal acceleration, or to control the average acceleration of a trajectory realization t = (y 1 , ..., y ℓ ). Again, the model can be adapted to both, but for practical reasons we will choose to control the maximal acceleration which involve simpler computations, thus leading to the measurement function Definition 3 (Maximal acceleration) The maximal acceleration of a tuple of points t = (y 1 , ...,

y ℓ ) is a max (t) = max 2 i ℓ-1 y i+1 -2y i + y i-1 .
(3)

We will now of course consider only trajectories having at least 3 points. Let T be a trajectory, Proposition 1 tells us that we can define a Number of False Alarms by an appropriate weighting of the probability P(a max (X T ) δ ). In the naive model, this probability only depends on the length ℓ of the trajectory (and, of course, δ ), and verifies Proposition 2 (Probability bound) If X T is a random trajectory of length ℓ, and δ is a positive real number, then

P(a max (X T ) δ ) (π • δ 2 ) ℓ-2 .
(4)

Proof -Assume that T = (k 0 , i 1 , ..., i ℓ ), and call B(x, r) the closed disc with center x and radius r.

Writing X ′ p = X k 0 +p-1 i p , we get P(a max ( X T ) δ ) = P ℓ p=3 X ′ p ∈ B(2X ′ p-1 -X ′ p-2 , δ ) ℓ ∏ p=3 P X ′ p ∈ B(2X ′ p-1 -X ′ p-2 , δ ) | X ′ p-1 , X ′ p-2 (π • δ 2 ) ℓ-2
because the area of B(x, δ ) ∩ [0, 1] 2 is bounded from above by π • δ 2 for all x. ✷ By Remark 1, we know that we can use the upper bound (4) to construct a NFA, as in (2). There are many possibilities to define the weights (w T ) T subject to the constraint ∑ T 1/w T 1. This gives us a way to adjust the detection thresholds for each structure. We choose to group trajectories together according to their length, dividing the set of trajectories into K groups T = T 1 ∪ ... ∪ T K (here, T ℓ denotes the set of trajectories of length ℓ), and weigh trajectories of a group uniformely by w T = K • #T ℓ for any T ∈ T ℓ . This choice is analyzed in detail in Section 2.6.3.

Proposition 3 (Continuous NFA for trajectories without holes) The family of functions (NFA T ) T ∈T defined by

∀ℓ, ∀T ∈ T ℓ , NFA T (δ ) = K(K -ℓ + 1)N ℓ • (π • δ 2 ) ℓ-2 (5)
is a Number of False Alarms for the measurement a max .

Proof -Since #T ℓ = (K -ℓ + 1)N ℓ , we have ∑ T ∈T 1 K(K -ℓ + 1)N ℓ = K ∑ ℓ=1 ∑ T ∈T ℓ 1 K • #T ℓ = 1,
and thus (5) defines a NFA thanks to Proposition 1. ✷

Let us quickly comment Proposition 3. We can rewrite (5) into NFA T (δ ) = K(K -ℓ + 1)N 2 • α ℓ-2 by using the relative density α = Nπδ 2 (which corresponds to the average number of points falling in a disc with radius δ ). We see that for a trajectory to be meaningful, we need to have α < 1. In other terms, only trajectories with maximal acceleration δ < 1/ √ Nπ might be detected as meaningful. Such kinds of bounds will be analyzed more precisely in Section 2.6.

Data quantization

In many applications, point detection is realized on a discrete grid of integer pixel coordinates, so that it may happen that three successive points in the sequence have a null acceleration. This is a very strong contradiction to the naive model, since this event has probability zero. Thus, if a long trajectory has a subtrajectory with a null acceleration, an algorithm that detects the most meaningful trajectories first will cut the longer trajectory into chunks to isolate the (optimal) null-NFA subtrajectory.

To avoid this kind of behavior we need to handle data quantization carefuly. There are at least two ways to do this: assume that the data have been properly quantized on the integer grid of the image and define a discrete version of the NFA, or consider a measurement imprecision and always consider the "worst-case scenario" for the measurements when computing accelerations.

First, we assume that the data have been quantized on an integer grid, say a rectangle Ω of Z 2 containing #Ω pixels. We can define a discrete version of the NFA by replacing the continuous acceleration area π •δ 2 by its discrete equivalent, the discrete acceleration area defined by

a d (δ ) = #S δ #Ω ,
where #S δ is the number of pixels enclosed in the discrete disc S δ = Z 2 ∩ B(0, |δ |) (see Fig. 2). In particular, when δ = 0 we obtain a discrete area of 1/#Ω that no longer leads to a null NFA. This leads to Definition 4 (Discrete maximal acceleration) The discrete maximal acceleration of a tuple of points t = (y 1 , ..., y ℓ ) is

a d max (t) = max 2 i ℓ-1 a d (y i+1 -2y i + y i-1 ). ( 6
)
Proposition 4 (Discrete NFA for trajectories without holes)

The family of functions

(NFA d T ) T ∈T defined by ∀ℓ, ∀T ∈ T ℓ , NFA d T (a) = K(K -ℓ + 1)N ℓ • a ℓ-2 (7)
is a Number of False Alarms for the measurement a d max .

Fig. 2 Discrete discs. A continuous disc and its corresponding discretization composed of all pixels whose centers lie inside the continuous disc. A discrete measure of area is better suited to the analysis of quantized data which might result in observing degenerate null-area continuous disc (the corresponding discrete disc has a non-null area).

We now examine the "worst-case scenario" acceleration. We assume that we have an estimate η > 0 of the measurements imprecision (corresponding roughly to the radius of one pixel in the previous example). We keep the same NFA than in the continuous case (Equation 5), but we replace the measurement function by

a w max (t) = max 2 i ℓ-1
a w (y i , y i+1 , y i+1 ), where

a w (x, y, z) = max dx,dy,dz∈ B(0,η) (x + dx) -2(y + dy) + (z + dz) .
One easily shows that a w (x, y, z) = x -2y + z + 4η, and therefore this amounts to keeping the same measurement function as in the continuous case and replacing the NFA (5) with

NFA w T (δ ) = K(K -ℓ + 1)N ℓ • (π • (δ + 4η) 2 ) ℓ-2 . ( 8
)
We can see that a null-acceleration trajectory will be counted as an acceleration of 4η, thus incurring a penalty to all accelerations. This is why in practice we assume that the point positions are quantized in a grid (or force this quantization if needed), and use the discrete NFA.

Algorithm

In this section we consider the discrete NFA given in Equation (7). When a meaningful trajectory is present, any slight deviation from it (removing or adding a point, for instance) will usually also be meaningful, so that we expect to detect a large number of overlapping meaningful trajectories. Hence, we choose to detect the trajectories greedily, by iterating the following process:

1. compute the most meaningful trajectory, that is, the one having the smallest NFA; 2. remove its points from the data.

To compute the most meaningful trajectory in a sequence of points, that is, K images I 1 , ..., I K , each containing N points, we use a dynamic programming strategy. Indeed, we compute for each point x in image I k the most meaningful trajectory ending in this point (note that in the following, we shall sometimes write x k instead of x to recall that x belongs to I k ). Denoting by G (x k , y k-1 , ℓ) the minimal acceleration of a trajectory of length ℓ ending with points y k-1 and x k , we obtain a Bellman equation

G (x k , y k-1 , ℓ) = 0 if ℓ = 2, min z∈I k-2 G (x, y, z, ℓ) otherwise, (9)
where

G (x, y, z, ℓ) = max a d (x -2y + z) , G (y, z, ℓ -1) . ( 10
)
This recursive formulation translates to Algorithm 1.

Algorithm: compute G input : f 1 , ..., f K the sets of {n k N} 1 k K points contained in each frame output: G for 2 k K do for x in f k do for y in f k-1 do G (x, y, 2) ← a d (0) for 3 ℓ k do G (x, y, ℓ) ← +∞ for z in f k-2 do a ← max(a d (x -2y + z), G (y, z, ℓ -1)) G (x, y, ℓ) ← min(a, G (x, y, ℓ)) end end end end end return G
Algorithm 1: Dynamic programming computation of G . We start by computing the values of G (x k , y k-1 , ℓ) for k = 2, then k = 3, ... each time reusing the results of the previous round according to the recursive formulation given by ( 9) and (10).

Now let us write

NFA d ℓ (a) = K(K -ℓ + 1)N ℓ • a ℓ-2 , (11)
so that NFA d T (a) = NFA d ℓ (a) for any trajectory T with length ℓ. Since the function a → NFA d ℓ (a) is monotone, the most meaningful trajectory with length ℓ is the one having the least maximal acceleration. Hence we can use Algorithm 1 to compute the smallest NFA of the point sequence (Algorithm 2). Moreover, if B(x k , y k-1 , ℓ) represents the most meaningful trajectory with length ℓ ending with points y → x (where t → x denotes the concatenation of trajectory t and point x), we can write

B(x k , y k-1 , ℓ) = y → x if ℓ = 2, B(y, ẑ, ℓ -1) → x otherwise, ( 12
)
where ẑ ∈ arg min z G (y, z, ℓ -1) (strictly speaking, the most meaningful trajectory might not be unique, so we have to choose one most meaningul trajectory arbitrarily; to simplify the description we shall nevertheless use the term "the most meaningful trajectory" throughout the paper). Finally, for each x k , the most meaningful trajectory ending in x k is B(x k , ŷk-1 , l) where

(ŷ k-1 , l) ∈ arg min (y k-1 ,ℓ) NFA d ℓ G (x k , y k-1 , ℓ) .
An algorithm similar to Algorithm 1 can hence be used to compute a trajectory having the smallest NFA (see Algorithm 2). In practice, we choose an arbitrary tuple (x, y, ℓ) such that there is an optimal trajectory of length ℓ ending on points y → x, and we extract a trajectory by backtracking, each time selecting the predecessor that minimizes locally the maximal acceleration among all predecessors that lead to an optimal trajectory. Note that we could also extract the trajectory having the minimal average acceleration (for example) among all optimal trajectories, simply by rerunning a similar dynamic-programming algorithm restricted to optimal trajectories. In the end, by applying the process greedily (as mentioned above), we obtain an algorithm that extracts all meaningful (or ε-meaningful) trajectories from a sequence of points (Algorithm 3).

It is often possible to save computation time by extracting several trajectories at once without recomputing the function G each time, because removing points from the current data set cannot decrease any value of G . This idea can be used, for instance, when the two most meaningful trajectories in the sequence do not share any point. The inner loop of Algorithm 4 does just that: if we have a way to define a set of disjoint minimal-NFA trajectories (for instance, greedily), we can extract them all at once (since they are of minimal NFA, and non-overlapping, we could extract them sequentially with interleaving G function recomputations, but this would not change the NFA of those trajectories, that would still be minimal).

Algorithm: minimal NFA input : G output: m the minimal NFA of a trajectory m ← +∞ for 3 k K do for x in f k do for y in f k-1 do for 3 ℓ k do m ← min(m, NFA d ℓ (G (x, y, ℓ))) end end end end return m
Algorithm 2: Find the minimal NFA value among all trajectories.

Algorithm: trajectory detection input : ε the maximal allowed NFA f 1 , ..., f K the sets of points contained in each frame output: S = t 1 , ...,t m the extracted trajectories

S ← / 0 repeat compute G m ← minimal NFA if m ε then t ← a trajectory of NFA = m S ← S ∪ {t},
and remove all points in t from the corresponding frames end until m > ε or there are no more points return S Algorithm 3: Greedy trajectories extraction using the NFA criterion.

Algorithm: trajectory detection accelerated input : ε the maximal allowed NFA f 1 , ..., f K the sets of points contained in each frame output: S = t 1 , ...,t m the extracted trajectories

S ← / 0 repeat compute G m ← minimal NFA stop ← false while m ε and stop = false do U ← { x, ∃ a traj. with NFA = m ending in point x } V ← a set of disjoint trajectories of NFA = m ending on a point of U S ← S ∪V
remove all points from the trajectories in V stop ← true if not all points of U have been removed m ← minimal NFA() end until m > ε or there are no more points return S Algorithm 4: Greedy trajectories extraction, accelerated by extracting several trajectories at once. Then if we have been able to extract a set of disjoint trajectories of minimal NFA that covers every pair of points where a trajectory of minimal NFA could end, we can continue doing the extraction for the next minimal NFA without recomputation. When this is no longer possible (because of some point removal) we need to recompute the G function to reactualize the NFAs.

We now examine the space and time complexity of algorithm 3 for an extraction round. The most expensive computation is that of function G . The space (memory) required is O(N 2 K 2 ), since we have to store a value for each triplet (x k , y k-1 , ℓ) in each image frame k. Each value computation takes O(N) operations because we have to consider all the points in the previous image, leading to a O(N 3 K 2 ) time complexity. The search for the minimal NFA and the extraction of the most meaningful trajectory have negligible time and space complexities. As the extraction must be repeated as long as there is any remaining meaningful trajectory, the global time complexity is O(sN 3 K 2 ), where s denotes the number of extracted ε-meaningful trajectories. In practice, on a standard PC desktop, for K = 50 images, the number N of points per image can reach several hundreds (and about one thousand for K = 20).

Variable number of points

In real data, the number of points is hardly ever constant throughout the sequence, so that instead of having N points in each of the K images, we have N 1 , N 2 , . . . N K points on images 1, 2, . . . K. Since the NFA is an upper bound on the average number of false alarms (Remark 1) we can simply take N = max k N k and keep the NFA unchanged. However, to obtain more accurate results, we can refine Proposition 4 with Proposition 5 (Discrete NFA for trajectories without holes, general case) The family of functions (NFA d T ) T ∈T defined, for any trajectory T with length ℓ starting at frame k 0 , by

NFA d T (a) = K(K -ℓ + 1) ∏ k 0 k k 0 +ℓ-1 N k • a ℓ-2 , ( 13
)
is a Number of False Alarms for the measurement a d max .

The proof is very similar to that of Proposition 3, except that now the set of trajectories T ℓ itself has to be decomposed with respect to the index of the starting frame.

Theoretical analysis

We now examine some theoretical consequences of Equation (5), which can have very practical consequences in the design of the data acquisition process. For simplicity reasons, we use the continuous NFA formulation, with a fixed number of points per image, but (7) or ( 13) would lead to similar conclusions.

Relation between the number of points and the maximal acceleration

Consider, as earlier, a sequence of K frames, each containing N points. We recall that the Number of False Alarms associated to a trajectory T with length ℓ 3 and maximal acceleration δ is (see Equation 5)

NFA T (δ ) = K(K -ℓ + 1)N ℓ • (π • δ 2 ) ℓ-2 .
Such trajectory is ε-meaningful as soon as

K(K -ℓ + 1)N ℓ • (π • δ 2 ) ℓ-2 ε,
which can be rewritten δ δ c , where the upper bound is the critical acceleration

δ c = 1 √ Nπ ε K(K -ℓ + 1)N 2 1 2ℓ-4 . (14)
Hence, as we already remarked at the end of Section 2.2, a necessary condition for trajectory detection is δ ∆ with ∆ = 1 √ Nπ , which gives an order of magnitude of the typical accelerations that can be handled by the NFA approach (and, in some sense, by any approach since accelerations greater than ∆ would allow detections in pure noise). Since the acceleration is inversely proportional to the squared frame rate (by doubling the frame rate, one divides accelerations by 4), this absolute bound can be useful in the design of the data acquisition process. Indeed, given the expected number of detected points in each frame (N), and the expected physical accelerations of objects (δ ), one can compute the critical frame rate, under which no trajectory detection is possible. Note, however, that the upper bound ∆ is not very accurate (see Table 1), thus using the exact value δ c (Equation 14) is probably a better idea. Table 1 Acceleration bounds ∆ = 1 √ Nπ and δ c (Equation 14) in function of N, expressed in pixel.image -2 in a 1000 × 1000 image for some values of l and K (ε = 1).

Influence of the trajectory length

A nice property of the a-contrario approach is that it permits us to relate different parameters by observing the way they are linked in the NFA formula. Table 1 shows that the trajectory length has a significant impact on the critical acceleration δ c (whereas the number of frames, K, has a much smaller impact). Thus, it could be interesting to study more precisely how the NFA balances the trajectory length and the acceleration, that is, how the critical acceleration grows as the trajectory lentgh increases. Since 1

K -ℓ + 1 K, we can write log(K(K -ℓ+1)) = 2β (ℓ) log K for some function β (ℓ) (taking values in [1/2, 1]), so that from (14) we get log δ c = log ∆ - log N + β (ℓ) log K -log √ ε ℓ -2 .
Hence, log δ c grows approximately like 1/ℓ, and attains for ℓ = K a value close to (and below) log ∆ . This is illustrated in Fig. 3. Last, we show the monotony of the critical acceleration with respect to the trajectory length (that is, the longer the trajectory, the looser the constraint on the acceleration).

Proposition 6 If ε 1, then the critical acceleration δ c given by ( 14) increases with respect to ℓ.

Proof -Rather than using ( 14), we go back to (5) and write, for ℓ ∈ {1, . . . , K},

F(ℓ, δ ) = N 2 K(K -ℓ + 1) • (Nπδ 2 ) ℓ-2
(thus, F(ℓ, δ ) is the NFA associated to a trajectory with size ℓ and maximal acceleration δ ). Now, if ℓ is such that F(ℓ, δ ) is smaller than 1 (since we assumed ε 1), then Nπδ 2 < 1 so that both ℓ → (Nπδ 2 ) ℓ-2 and ℓ → N 2 K(K -ℓ + 1) are decreasing with respect to ℓ, and so is ℓ

→ F(ℓ, δ ). Hence, if ℓ 1 > ℓ 2 , we have, for δ = δ c (ℓ 2 ), F(ℓ 1 , δ c (ℓ 2 )) < F(ℓ 2 , δ c (ℓ 2 )) = ε, which proves that δ c (ℓ 1 ) > δ c (ℓ 2 ) since δ → F(ℓ 1 , δ ) is in- creasing.
✷

Asymptotic bounds and the importance of the combinatorial factor

Now we would like to assess the importance of the combinatorial factor in the definition of the NFA. As was discussed above before Proposition 3, there are several ways to define the weights of the structures. In (5), we chose to weigh the trajectories uniformely with respect to their length (that is, such that the expected number of false alarms is equally shared among all possible trajectory lengths), that is

NFA T (δ ) = K(K -ℓ + 1)N ℓ • P(a max (X T ) δ ). ( 15
)
Another more classical choice would be to set a uniform weight w T = #T for all trajectories, thus obtaining

NFA ′ T (δ ) = K ∑ m=1 (K -m + 1)N m • P(a max (X T ) δ ). (16)
Suppose that we observe a trajectory t with length ℓ and maximal acceleration δ = a max (t), what difference will each NFA definition make? This trajectory is detected if NFA T (δ )

(or NFA ′ T (δ )) is below a certain threshold ε, hence it is in- teresting to estimate the ratio NFA ′ T (δ ) NFA T (δ ) = 1 K(K -ℓ + 1)N ℓ K ∑ m=1 (K -m + 1)N m 1 K 2 N ℓ K ∑ m=ℓ N m N K-ℓ K 2 .
This lower bound shows that when ℓ is small, the detection penalty incurred when using NFA ′ is very large and it will thus be more difficult to detect small trajectories with this criterion.

In the following, we compare more precisely NFA and NFA ′ and compute asymptotic estimates when the number of frames (K) becomes large. Let us deal first with the function NFA T (δ ). We consider a trajectory T spanning #T = µK images of the sequence for a fixed µ ∈ (0, 1], with a maximal acceleration δ , among N points per frame. We write α = Nπδ 2 , and propose to determine, when K gets very large, if the trajectory is meaningful. We have

NFA T (δ ) = K • (K -µK + 1) • N µK • (πδ 2 ) µK-2 ∼ K→∞ π -2 δ -4 K 2 (1 -µ)α µK , so that log NFA T (δ ) ∼ K→∞ µK log α and lim K→+∞ NFA T (δ ) 1 ⇐⇒ α < 1 (17)
This means that for any µ ∈ (0, 1) and any maximal acceleration δ < ∆ (that is, such that α < 1), for K large enough, any trajectory with maximal acceleration at most δ that spans µK frames among K will be meaningful. In practice, trajectories that lead to values of α near 1 are difficult to detect because they need to be observed on a large number of images (large K) to be meaningful. This phenomenon is illustrated in Fig. 4. Another non-asymptotic relationship between the acceleration δ , the number of frames K and the number of points per frame N is illustrated in Fig. 5. Fig. 4 Non-asymptotic counterpart of (17). Theses three curves represent, as a function of α = Nπδ 2 , the minimum number of frames (K) required for a trajectory with length ℓ = ⌊µK⌋ and maximal acceleration δ to be 1-meaningful according to the NFA criterion (15). Upper (blue) curve: µ = 0.3; middle (green) curve: µ = 0.5; lower (red) curve: µ = 1. The number of points per frame is N = 100. Now we study the asymptotic behavior of NFA ′ T (δ ) (see Equation 16). First, we notice that

K ∑ m=1 (K -m + 1)N m = N K K ∑ m=1 mN -m-1 ∼ K→∞ N K (N -1) 2 and thus NFA ′ T (δ ) ∼ K→∞ N K+2-µK α µK-2 (N -1) 2 . Therefore, log NFA ′ T (δ ) ∼ K→∞ K (1 -µ) log N + µ log α
(δ ) < 1 ⇐⇒ log α < µ -1 µ log N. ( 18
)
Since µ-1 µ → -∞ as µ → 0 + , Equation ( 18) shows that it is indeed asymptotically much harder to detect small trajectories with NFA ′ than with NFA. This fact is illustrated by Fig. 6, on which we can see that even when µ is slightly smaller than 1, NFA permits to detect more trajectories than NFA ′ , both asymptotically and non-asymptotically. Each curve represents, as a function of µ, the maximal value of α = Nπδ 2 allowed for a trajectory with length ℓ = ⌊µK⌋ and maximal acceleration δ to be 1-meaningful in a sequence of K images, among N = 100 points per frame. The red upper curves are obtained with the criterion NFA, whereas the blue lower curves are obtained with the criterion NFA ′ . The full curves correspond to K = 100, and the dashed curves correspond to the asymptotical estimates obtained when K → +∞, that is, log α = 0 for NFA and log α = log N • (µ -1)/µ for NFA ′ . These curves clearly demonstrate that not only NFA is better suited to the detection of small trajectories than NFA ′ (it allows for trajectories having a much larger maximal acceleration), but it is also more efficient even for relatively large values of µ (NFA ′ being slightly better only for almost complete trajectories). Asymptotically, NFA is always the best choice.

Trajectories with holes

In many practical situations, because of occlusions or acquisition noise, some trajectory points will not be detected in one or more frames. In this section, we generalize the previous framework to trajectory detection in the case of missing points.

Number of false alarms

The naive model remains unchanged (keeping the notations of the previous section): we are given K images I 1 , ..., I K , each image I k containing N points X k 1 , ..., X k N , and we assume by Helmholtz principle that all random points X k i are independent and uniformly distributed in [0, 1] 2 .

Definition 5 (Trajectories with holes) A trajectory of size s is a sequence of pairs T = {(i 1 , τ 1 ), ..., (i s , τ s )}, such that τ 1 < ... < τ s . We denote by T the set of all trajectories, and by T s the set of trajectories of size s. We bijectively associate to the trajectory T the tuple of random (i.i.d., uniformely distributed) variables X T = (X τ 1 i 1 , ..., X τ s i s ), that we also abusively call a (random) trajectory.

Since the definition above is more general than the special case of trajectories without holes (Definition 2), we chose to keep the same word (trajectory) and the same notations (T, T s ) as in Section 2. This should not lead to ambiguities, since we will only consider trajectory with holes in this section.

As in Section 2, we would like to build an a-contrario detection model to detect trajectories (here, with holes). We consider three parameters of interest for the computation of the NFA of a trajectory X T = (X τ 1 i 1 , ..., X τ s i s ): the trajectory length, its size and its number of runs. The length is the total number of frames that the trajectory spans (τ s -τ 1 + 1), the size is the number of (detected) points it contains (s), and a run is a maximal set of consecutive points. Note that if we call hole a maximal set of consecutive missing points, then the number of runs equals the number of holes plus one.

We first need to generalize the notion of maximal acceleration a max (Definition 3) to the case of trajectories with holes. A natural way to do this consists in interpolating the missing points of the trajectory and compute its maximal acceleration. Since we would like to keep using an algorithm based on dynamic programming, we use the most local choice, that is, a simple constant speed interpolation. This leads to the following Definition 6 (Maximal acceleration with holes) The maximal acceleration of the realization t = (y τ 1 1 , ..., y τ s s ) of a trajectory T is

a h max (t) = max 2 i s-1 a h (y τ i-1 i-1 , y τ i i , y τ i+1 i+1 ) , (19)
with, for all points x i , y j , z k (i

< j < k), a h (x i , y j , z k ) = z -y k -j - y -x j -i . ( 20
)
We now compute, as in Proposition 2, a probability bound for the maximal acceleration of a random trajectory with holes.

Proposition 7 (Simple probability bound) If a random trajectory X T with size s has holes of size h 1 , ..., h p-1 , then for any δ > 0 one has

P(a h max (X T ) δ ) (π • δ 2 ) s-2 • ∏ 1 i p-1 (h i + 1) 2 . ( 21
)
Proof -We assume that T = {(i 1 , τ 1 ), ..., (i s , τ s )}, and write

X ′ q = X τ q i q and M q = X ′ q-1 + τ q -τ q-1 τ q-1 -τ q-2 (X ′ q-1 -X ′ q-2 ), so that P(a max ( X T ) δ ) P s q=3 X ′ q ∈ B(M q , (τ q -τ q-1 )δ ) s ∏ q=3 P X ′ q ∈ B(M q , (τ q -τ q-1 )δ ) | X ′ q-1 , X ′ q-2 (π • δ 2 ) s-2 s ∏ q=3 (τ q -τ q-1 ) 2 (π • δ 2 ) s-2 p-1 ∏ i=1 (h i + 1) 2 .

✷

For efficiency reasons, we want to design an algorithm that can share computations, that is, we want to be able to reuse the computations made on subtrajectories and extend them to obtain the results for bigger trajectories. To do this efficiently, we shall not consider, for a given trajectory, the individual sizes h i of its holes, but simply its length ℓ, its size s and its number of runs p. This is why we derive from (21) the following Proposition 8 (Practical probability bound) If a random trajectory X T has length ℓ, size s and number of runs p, then for any δ > 0 one has

P(a h max (X T ) δ ) (π • δ 2 ) s-2 • ℓ -s p -1 + 1 2p-2 (22)
with the convention that the right-hand parenthesis equals 1 (( 0 0 + 1) 0 ) when p = 1.

Proof -We consider the maximum value of the right-hand term of ( 21) over all possible hole sizes h 1 , ..., h p-1 that are feasible for parameters ℓ, s and p. Relaxing the constraint that the h i have to be integers, we face the optimization problem max

(h i ) p-1 ∏ i=1
(h i + 1) 2 ; ∑ i h i = ℓs, and ∀i, h i 0, which, denoting ξ i = h i + 1, has the same solutions as the problem max

ξ ∈C E(ξ ), with E(ξ ) = p-1 ∑ i=1 log(ξ i ) and C = ξ ∈ [1, +∞) p-1 , p-1 ∑ i=1 ξ i = ℓ -s + p -1 .
Now if ξ ∈ C has not identical coordinates, we can choose two different values, say 1 ξ i 1 < ξ i 2 , and replace them both by

(ξ i 1 + ξ i 2 )/2 1 to form a new ξ ′ ∈ C. Then, E(ξ ′ ) -E(ξ ) = 2 log((ξ i 1 + ξ i 2 )/2) -log(ξ i 1 ) -log(ξ i 2 ),
and this quantity is positive by strict concavity of the log function. Thus, the unique solution of max ξ ∈C E(ξ ) satisfies ξ i = (ℓs)/(p -1) + 1 for all i, and the maximum value of ∏ i (h i + 1) 2 over feasible hole sizes h 1 , ..., h p-1 is bounded from above by

ℓ -s p -1 + 1 2p-2
.

Using this bound in (21) yields the announced result. ✷

Now we need to define the combinatorial term w T that premultiplies the NFA. Recalling the discussion of the end of Section 2.2, we choose to group the trajectories by their length and size, and to use uniform weights in each category. The number of trajectories of length ℓ and size s that can fit in K frames is bounded from above by (K -ℓ + 1) ℓ s N s , and since we cluster the trajectories by their lengths and sizes, we have to count the number of such clusters. Indeed, it is bounded from above by Kℓ, since there are less than K ways to choose the length, and knowing that the length is ℓ there are less than ℓ ways to choose the size. Combining these remarks with Proposition 8 establishes the following Proposition 9 (NFA for trajectories with holes) The family of functions (NFA T ) T ∈T defined for any trajectory T of length ℓ, size s and number of runs p by

NFA T (δ ) = Kℓ(K -ℓ+1) ℓ s N s (πδ 2 ) s-2 ℓ -s p -1 + 1 2p-2 (23)
is a Number of False Alarms for the measurement a h max .

This new function NFA T is a kind of generalization of (5). Indeed, for a trajectory T without hole (that is, such that p = 1, and consequently ℓ = s), we have

NFA T (δ ) = ℓ • K(K -ℓ + 1)N ℓ (πδ 2 ) ℓ-2
which is, up the a new factor ℓ, the value given in (5). This new factor simply comes from the fact that we do not know a priori that the number of runs of the trajectory is one.

Algorithm

In the practical implementation that we describe below, we use a discrete version of the acceleration, obtained by replacing the norm involved in ( 19) by a discrete area measure, exactly as we did in Section 2 (see Definition 4).

We want to compute, for each point x of each image I i (that we denote by x i ), the most meaningful trajectory that ends in x i (or, to be more precise, one of such most meaningful trajectories) . This information can be extracted from the function G (x i , y j , ℓ, s, p), which represents the least maximal acceleration of a trajectory of length ℓ, size s, and having p consecutive runs (that is, p -1 holes), ending with the point y j in frame j < i followed by the point x i in frame i.

We say that a tuple (i, j, ℓ, s, p) is undefined if there is no trajectory ending with its two last points in frame i and j, with length ℓ, size s and having p runs of consecutive points. For instance, if ℓ < s or s < 2, the tuple is undefined. We define for i > j

G (x i , y j , ℓ, s, p) =    +∞ if (i, j, ℓ, s, p) is undefined, 0 if s = 2, min u 1,z∈I j-u Ḡ (x, y, z, ℓ, s, p) else, ( 24
)
with the convention, for i

> j > k, that Ḡ (x i , y j , z k , ℓ, s, p) = max a h (x, y, z), G (y, z, ℓ -(i -j), s -1, p -1 i = j+1 ) (25)
and as usual 1 a =b = 1 if a = b and 0 otherwise. Notice that as in Section 2.4, the superscript k in x k simply reminds us that the point x belongs to image I k , so we sometimes omit it and simply write x.

We deduce from (25) a dynamic programming algorithm to compute G , similar to the one we presented in Section 2 for the trajectories without holes. We can then backtrack to find the most meaningful trajectory ending in each point x i by defining, for i > j, the recursive function

B(x i , y j , ℓ, s, p) =    undefined if (i, j, ℓ, s, p) is undefined, y → x if s = 2, B(y, ẑ j-û, ℓ -(i -j), s -1, p -1 i = j+1 ) → x else,
where ẑ j-û realizes the minimum in the last line of (24). Finally, for each x i , the most meaningful trajectory with holes ending in x i is B(x i , ŷ ĵ, l, ŝ, p), where ( ĵ, ŷ ĵ, l, ŝ, p) = arg min j<i,y j ,ℓ,s,p NFA d ℓ,s,p (G (x i , y j , ℓ, s, p)).

We can analyze the spatial and temporal complexities of the algorithm. As in the case of trajectories without holes, the most costly operation is the computation of the G function. Its complexity is O(N 2 K 5 ) in space and O(N 3 K 6 ) in time. Since the extraction must be repeated until there are no more meaningful trajectories, the global time complexity is O(sN 3 K 6 ), where s is the number of extracted ε-meaningful trajectories. In practice, on a standard PC desktop, for K = 30 images, the number N of points per image can go up to about one hundred. For long sequences containing much more than 30 images, the algorithm cannot be used directly (due to the N 6 term in the time complexity), but one could probably obtain good results by cutting the image sequence into small parts and applying the algorithm on each part (raising the issue of long trajectories spanning several parts).

Variable number of points and rectangular images

As for the NFA of trajectories without holes (Section 2), we can adapt the NFA given in Proposition 9 to the case of a variable number of points per image. Let us write N k the number of points present in image k. The simplest strategy consists in applying directly the definition of Proposition 9 with N = max k N K . If N k has strong variations, a more sensitive detection can be obtained by replacing in NFA T (δ ) the term N s by max

k 0 =i 1 <i 2 <...<i s =k 0 +ℓ-1 N i 1 • ... • N i s ,
where T is a trajectory starting in image k 0 , with length ℓ and size s. This term is easily computed once the sequence

(N k ) 1 k K has been sorted.
As in the case of trajectories without holes, the NFA can also be adapted to rectangular images (see Section 2.3).

Theoretical results

We now analyze the asymptotic behavior of the NFA on some particular cases. They are all composed of one trajectory spanning the K images, and N -1 additional spurious points in each frame. The trajectory has a maximal acceleration of δ .

Long trajectory with a single hole

We first study the case where the trajectory is composed of two parts separated by a unique hole of length h = εK. We thus have ℓ = K, s = (1 -ε)K, p = 2, and we write α = Nπδ 2 as in Section 2.6.3.

We study under which conditions, when K gets large, the trajectory is meaningful, and if it is more meaningful than its first (or equivalently last) part. First we derive an asymptotic expansion of

NFA T (δ ) = K • (K -ℓ + 1) • ℓ • ℓ s • α s (πδ 2 ) 2 • (εK) 2(p-1) = K 2 K (1 -ε)K α (1-ε)K (πδ 2 ) 2 (εK) 2 .
From Stirling's Formula, one easily derives the expansion

log K ηK = -Kh(η) - 1 2 log K + O K→+∞ (1),
where η ∈ (0, 1) is fixed, and

h(η) = η log(η) + (1 -η) log(1 -η).
Hence, we have

log NFA T (δ ) = K (1 -ε) log α -h(ε) + 7 2 log K + O K→+∞ (1),
which proves the asymptotic equivalence lim

K→+∞ NFA T (δ ) 1 ⇐⇒ log(α) < h(ε) 1 -ε . ( 26
)
This asymptotic condition on α is illustrated in Fig. 7. We can notice that the asymptotic limit on log α given by the righ-hand term of ( 26) is quite accurate (and almost linear) as long as the relative hole size ε is not too large. If the hole size is half the trajectory length (ε = 1/2), then the asymptotic condition is α < 1 4 . Now we would like to investigate the condition under which the complete trajectory is more meaningful than its starting or ending parts. Writing γ = (1 -ε)/2 so that each small trajectory has a size γK, this condition writes

NFA T 1+2 (δ )/NFA T 1 (δ ) 1, that is K K • 1 (1 -γ)K + 1 • K γK • K 2γK γK γK • α 2γK α γK • ((1 -2γ)K + 1) 2 1 1 or equivalently -Kh(2γ) + γK log α + log K + O K→+∞ (1) 0.
hole relative length (ε) 0 0.2 0.4 0.6 0.8 1 1 10 -2 10 -4 α Fig. 7 Influence of the hole size. Plot, in function of ε, of the maximal value of α (in log scale) for a trajectory with a single hole of size h = ⌊εK⌋ to be meaningful. The total number of points in each frame is N = 100, and the length of the sequence is K = 100 for the green lower curve and K = 400 for the blue middle curve. The red upper curve is the asymptotical limit h(ε)/(1 -ε) of log α corresponding to the case K = +∞. We can see that when the trajectory hole becomes fairly important, the maximal allowed acceleration for a trajectory to be meaningful plummets.

Since h(2γ)/γ = 2h(ε)/(1 -ε), we thus have the property that the whole trajectory is asymptotically more meaningful than its parts when K → +∞ if and only if

log(α) < 2 h(ε) 1 -ε . (27)
Note that this inequality constraint on α is stronger (the quantities are negative) than the one obtained in (26), hence the case when a trajectory is 1-meaningful but less meaningful than its parts can be encountered.

Dotted trajectories

If the trajectory is made of a succession of single points and one-frame holes, what is the condition as K → +∞ to have a meaningful trajectory? We now have ℓ = K, s = (K + 1)/2 (K being odd), p = (K -1)/2, so that from (23) we can derive the asymptotic expansion

log NFA T (δ ) = K 2 log(16α) + 3 2 log K + O K→+∞ (1). (28)
Hence, a dotted trajectory is asymptotically meaningful when K → +∞ as soon as α < 1 16 .

(29)

Dashed trajectories

In the more general setting of a dashed trajectory made of p runs of u consecutive points separated by holes spanning v frames, we have ℓ = K = p(u + v)v and s = pu, so that if u and v are fixed,

log NFA T (δ ) = 2 log p + log p(u + v) -v pu + pu log α + 2p log(1 + v) + O p→+∞ (1). Now we have pu p(u + v) -v = η + O p→+∞ 1 p with η = u u + v , so that log p(u + v) -v pu = -p(u + v)h(η) - 1 2 log p + O p→+∞ (1)
and

log NFA T (δ ) = 3 2 log p+ pu log α + 2 u log(1 + v) - h(η) η + O p→+∞ (1).
Hence, a dashed trajectory is asymptotically meaningful when p tends to infinity if and only if

log α h(η) η - 2 u log(1 + v), (30)
where η = u u+v is the asymptotic density of known points. This formula yields an interesting relation between the density of known points and the allowed maximum acceleration, as illustrated in Fig. 8. When u = v = 1, we have η = 1/2, h(η) =log 2 and the right-hand term of (30) islog 16, in accordance with (29).

The ROADS algorithm, and comparison

The ROADS tracking algorithm

As we mentioned in Introduction, the state-of-the-art ROADS algorithm (Veenman et al, 2003b) is an extension of GOA that compares favorably to most classical tracking algorithms like MHT [START_REF] Reid | An algorithm for tracking multiple targets[END_REF], [START_REF] Salari | Feature point correspondence in the presence of occlusion[END_REF], Rangarajan and Shah (1991), [START_REF] Chetverikov | Feature point tracking for incomplete trajectories[END_REF], and GOA itself [START_REF] Veenman | Resolving motion correspondence for densely moving points[END_REF]. It can handle points entering and leaving the scene, as well as missing and spurious points. It requires the setup of several parameters, which are listed in Table 2.

The criterion measuring the local smoothness of a trajectory on the consecutive points (x, y, z) is

φ (x, y, z) = w 1 - v(x, y) • v(y, z) v(x, y) • v(y, z) + (1 -w) 1 -2 v(x, y) • v(y, z) v(x, y) + v(y, z)
density of points (η) 0 0.2 0.4 0.6 0.8 1 1 10 -2 10 -4 10 -6 α density of points (η) 0 0.2 0.4 0.6 0.8 1 1 10 -2 10 -4 10 -6 α Fig. 8 Required precision for "dashed" trajectories. Left: The green lower curve shows the maximal value of α allowed (in log scale) for a trajectory to be meaningful when it is made of repetitions of runs of length u and holes of length v, with u + v = 20 and u = ⌊η(u + v)⌋. The length of the sequence is K = 100, and the total number of points in each frame is N = 100. The upper red curve is the asymptotical limit h(η)/η -2 log(1 + v)/u corresponding to p = +∞. (note that the staircasing effect is due to the definition of u). As long as the density of known points (η) is large enough, the critical value of log α is quite well approximated by the asymptotic bound and the relation to η is almost linear. When the density η becomes too small, the maximal allowed acceleration for the trajectory to be meaningful quickly plummets.

Right: The asymptotic condition on α given by ( 30), for a variety of values of u+v = 10 (lower green curve), 20 (middle blue curve) and 40 (upper red curve). The curves are close, hence showing that the minimal required precision for dashed trajectories to be meaningful mostly depend on the density η = u/(u + v) and not much on the period u + v of the runs. where v(x, y) = yx. As we can see, this criterion combines (with a weight parameter w) an angular variation (first term) and a speed variation (second term). Assume that M k objects are tracked until the k th frame, and N k+1 points are observed in frame k +1. The trajectories already constructed can either link to one of the observed points, or to a missing "slave measurement", meaning the corresponding object in frame k + 1 is missing. Additionally, a point of frame k + 1 can be tagged as spurious. All these possibilities are called individual assignments.

w
Each individual assignment a has a cost c(a). The cost of linking a trajectory to a point in frame k +1 is the smoothness criterion as defined above (if one of the past measurements is missing, we estimate its position through linear interpolation). The cost of considering a point in frame k + 1 as spurious and the cost of linking a trajectory to a slave (missing) measurement are both equal to the value of the parameter φ max (eg. a missing point has the cost of the worst possible trajectory continuation). The algorithm restricts its possible assignments using its cut-off values, for instance, two points in consecutive frames can be linked only if they are at most d max pixels appart.

The local cost of all the individual assignments between two consecutive frames is obtained by averaging their costs. Let A k = {a 1 , ..., a p } be the set of individual assignments between frame k and k + 1, that is, such that every trajectory tracked in frame k appears in exactly one of the assignments, and every measurement in frame k +1 appears in exactly one of the assignments, then the local cost is

C k (M k , A k ) = 1 M k p ∑ i=1 c(a i ).
The optimization of this cost for a fixed k is a minimumweight perfect matching problem, and can be solved efficiently using for instance the Hungarian algorithm [START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF]. Finally, the global motion model averages costs over the whole sequence, leading to the minimization of

C(A) = K-1 ∑ k=2 C k (M k , A k )
where A = (A 2 , ..., A K-1 ) is a multi-assignment. Other optimization objectives are: as many points as possible should be included in a trajectory, and there should be as few trajectories as possible. In its generality, the global motion model optimization is a NP-hard problem, thus intractable in practice. One of the approximation made by the ROADS algorithm is to sequentially optimize the global model on a restrained time window (typically using s = 2 or s = 3), hence computing

A k:s min = argmin A k:s C k:s (A k:s )
where

C k:s (A k:s ) = s ∑ p=1 C k+p-1 (M k+p-1 , A k:s [p])
and A k:s = (A k , ..., A k+s-1 ) is a multi-assignment. The approximation to the global solution is then

A = (A 2:s min [1], ..., A K-1:s min [1]).
This approach results in an initialization problem at the beginning of the sequence: the assignment between the first two frames is considered given. To mitigate this strong requirement, the ROADS algorithm uses a "minimal-motion" criterion c(x, y) = yx to initialize the assignment between the first two frames of the sequence, and then a successive up-and down-processing to reduce the imprecision of the initial assignments. We refer the reader to Veenman et al (2003b) for a detailed explanation, including how tracks can be ended and started.

The core of the ROADS algorithm (see Algorithm 5) is the computation of the local scope optimization. In order to compute A k:s at each frame k, the algorithm recursively enumerates all potential assignments between successive frames in the scope. Of course, this set of assignments is too large to be exhaustively enumerated, and therefore the algorithm uses a branch-and-bound approximation strategy. It first makes an "optimal cost bound" guess C b by initializing the global solution on the time scope with the local solutions of the minimum-weight perfect matching between each two consecutive frames in the scope. This cost bound is then gradually lowered.

At each recursion step (that is, each frame in the scope), the bound on the current optimal matching cost is lowered by using a cost-bound constraint called γ max . It is derived from the cost C k min = C k (M k , A k min ) of the best possible assignment A k min between frames k and k + 1 (which can be obtained by the Hungarian algorithm) and the current global cost bound C b by

γ max = min(F γ l C k min , F γ g C b /s ′ ),
where

F γ l 1 is the local cost factor, F γ g
1 is the global cost factor and s ′ is the length of the remaining scope.

The intuition behind this bound is that the cost of the assignment at frame k corresponding to the optimal solution on the time scope cannot be too far from the cost C k min of the optimal (local) assignment between the frames k and k + 1, and that the cost of the assignment on the time scope is more or less uniformly distributed between all pairs of frames, and thus should not be too far from C b /s ′ . To enumerate the successive best assignment between frames, ROADS uses Murty's algorithm [START_REF] Murty | An algorithm for ranking all the assignments in order of increasing cost[END_REF], which takes a cost matrix D k and a set of previous assignments Y and returns the next best assignment not in Y . The set of all assignements between frames k and k + 1 is denoted by U k .

The costMatrix(A k-1 , k) function returns a matrix containing the cost of each possible assignment between a trajectory of A k-1 and a point of the kth image.

Experiments

In the following, we propose to compare the NFA algorithm with ROADS to evaluate its strengths and weaknesses against a state-of-the-art solution. We start with experiments on synthetic data, similar to those used by the authors of ROADS in their presentation papers (Veenman et al, 2003a,b). Let us first briefly present the way they generate trajectories using the Point-Set Motion Generation (PSMG) algorithm (experiments having different parameters will be signaled):

Algorithm: ROAD(A k-1 , k,C b , A k:s sol ) input : A k-1 the previous assignment,
k the current frame number, C b the current cost bound, A k:s sol the current best assignment output: A k:s sol the new best assignment

D k ← costMatrix(A k-1 , k) if s = 1 then A k min = minCostAssignment(D k ) if C k (M k , A k min ) < C b then A k:s sol ← (A k min ) end else Y ← / 0 repeat A ← nextBestAssignment(Y, D k ) Y ← Y ∪ {A} C 0 ← C k (M k , A) T sol ← A k:s sol [2..s] R ← ROAD(A, k + 1, s -1,C b -C 0 , T sol ) A k:s = (A) :: R if C k:s (M k , A k:s ) < C b then C b ← C k:s (M k , A k:s ) A k:s sol ← A k:s end γ max = min(F γ l C k min , F γ g C b /s) until Y = U k or C 0 C b or C 0 γ max end return A k:s sol
Algorithm 5: Core of the ROADS algorithm the initial position of each trajectory is chosen uniformly at random in the first image; the initial velocity magnitude is chosen using a normal random variable v 0 ∼ N (µ = 5, σ = 0.5) and its angle β 0 is chosen using a uniform distribution in [0, 2π]; the velocity magnitude and angle are updated in each frame using

v k+1 ∼ N (µ = v k , σ = 0.2) β k+1 ∼ N (µ = β k , σ = 0.2).
The image domain is divided in 100 × 100 pixels, and the length of the sequence is set to 20 frames (see Fig. 9 for an illustration of the trajectories generated). Most of the experiments are realized with 20 trajectories (like in the ROADS paper).

Since the ROADS authors were comparing their algorithm with an algorithm that did not allow trajectories entering or leaving the scene, they required that all trajectories fit completely inside the frames and span the whole sequence, and we will usually do the same (if a trajectory does not fit inside the frame, we regenerate it). They also impose that in the experiments where points are missing, all points are still detected in the first and last two frames. To have experiments coherent with theirs we generally impose the same constraints, but in some experiments (with a great number of trajectories in the images) constraining trajectories to stay inside the frame seemed unnatural since it forced trajectories to have a beating and swirling motion inside the frames. We chose to keep a constant number of trajectories, but to allow them to leave the image, and when this happens, to generate a new random trajectory that starts on the border of the frames (we force all trajectories to have at least three points in the frame).

Fig. 9 Point-Set Motion Generation (PSMG) algorithm. Here we display a sample of 20 trajectories spanning 20 frames generated with the PSMG algorithm, which will be used to produce synthetic data to assess the performance of the ROADS and NFA algorithms. The trajectories have an homogeneous and smooth motion. The points of some trajectories have been highlighted to show the speed.

Additionaly, we impose the following constraints:

we quantize the trajectory coordinates to the nearest integers, thus implicitely defining the scale of the measurements to the size of one pixel; to avoid ambiguities when comparing the detection results to the ground truth, trajectories cannot share points (otherwise we regenerate one of the trajectories) and when adding noise points we avoid covering already existing points.

when we remove points we target solely trajectory points (we do not remove noise points, so we can make experiments with a varying number of trajectory points removed, while keeping a constant number of noise points). We choose a certain uniform probability α of removing a point.

Performance estimates gathered in the experiments below are averages of a measure over 400 runs of the algorithmson random data. However, in some experiment results, the measure that we compute might be undefined (for instance when no trajectory was detected). In this case, we only take experiment results that have a defined measure into account, and measurements might thus consist of averages of less than 400 repetitions.

A well-known interest of a-contrario methods are their small number of parameters, which simplifies their use and their study. More accurately, the NFA algorithm has exactly one explicit parameter, the maximal NFA value of a trajectory we can extract. The effect of this threshold is simple: it drives the selection of a subsequence of the successively extracted trajectories. In other words, if ε < η, T(ε) ⊆ T(η), where T(x) is the set of trajectories extracted by the algorithm for a maximal value of the NFA equal to x. This implies that changing the threshold will not dramatically change the results, contrary to methods like ROADS that use their parameters in the computations. In practice, as usually done in a-contrario methods (see [START_REF] Desolneux | From Gestalt Theory to Image Analysis: A Probabilistic Approach[END_REF]) and unless otherwise specified, we set this threshold ε to 1.

In contrast, ROADS has many parameters, which can be tuned to set ROADS in different "modes" that may be better suited to certain types of data. Since these parameters might (at least in theory) be learned on data, we felt it was fair to try several sets of parameters and show the best results that could be achieved in the comparisons.

Here is the way we proceeded: we tested six "modes" for the ROADS algorithm on a small batch of data (40 repetitions) for each experiment. We then selected the three modes that would compare the best with the NFA algorithm on the various experiments. Some of the modes will have strengths and weaknesses compared to others, but they mostly have the same global behaviour. In practice, the strengths and weaknesses of the NFA method when compared to ROADS do not dramatically change when including several modes, rather than just the most general parameters for the ROADS algorithm (mode 1 below). However, we include the results of the three selected modes for the sake of completeness.

To be fair with the ROADS algorithm, which relies on knowing the maximal speed and maximal smoothness criterion of the trajectories in the data, we compute these values and give them to the algorithm. More precisely, for each experiment, and each parameter (eg. number of noise points added), we compute the maximal speed d max and maximal smoothness criterion φ max before crippling (eg. removing points) across the batch of 400 repetitions (rather than on a per-file basis), and we feed them to the algorithm when processing those 400 repetitions.

The first mode is the general ROADS algorithm with the minimal number of present points set to p min = 1 and the maximal number of interpolated points equal to a max = +∞. The second mode is p min = 1 and a max = 0, that is, we disallow interpolation. The third mode is p min = 3 and a max = 0, which means that we disallow interpolation and we expect to see at least 3 consecutive points on each trajectory segment. For the three other modes, we set p min = 3, a max = 3, but rather than choosing the maximal speed and maximal smoothness criterion as given by their maximal value on the batch of 400 repetitions, we select in turn: d 4 max = 0.8 • d max , φ 4 max = 0.8 • φ max for mode 4, d 5 max = 0.5 • d max , φ 5 max = 0.5 • φ max for mode 5 and d 6 max = 0.5 • d max , φ 6 max = 0.8 • φ max for mode 6.

The other default ROADS parameters given in the implementation that was sent to us by its authors were kept unchanged (w = 0.1, F ℓ = F g = 1.05, s = 2). We tried to make some experiments with s = 3, but this would generally not change the results (and sometimes even have a negative impact) and be much more computationally intensive. The maximal number of hypotheses h max that ROADS can explore when trying to find the best assignment in the time scope has been kept equal to 300 as in the given implementation.

After running all the ROADS modes on small batches of repetitions, we selected modes 1 (original ROADS), 3 and 4 as giving the best results. They will now be called modes A, B and C (see Table 3). Note that these modes are not real algorithms, since their parameters depend on true data values (d max and φ max ) that are not estimated but computed from an oracle. In that sense, the methodology we use to compare the NFA and ROADS algorithms minimizes the issue of parameter selection that is recurrent with ROADS (but this issue will be discussed later, in particular in Section 4.3.4). Note also that ROADS results could probably be slightly improved by trying a larger number of parameters, however our goal is not to make a study of ROADS, but rather to give an idea of the state of the art performances, to make it possible for the reader to appreciate the NFA results.

mode p min a max d max φ max A 1 +∞ 1 • d max 1 • φ max B 3 0 1 • d max 1 • φ max C 3 3 0.8 • d max 0.8 • φ max
Table 3 Parameters defining the three (best) modes of the ROADS algorithm in the experiments used for comparison with the NFA algorithm. Note that these three modes are based on an oracle that observes the values of φ max and d max on the (supposedly unknown) true trajec- tories.

Comparison criteria

In the literature, tracking algorithms are generally compared using two sets of criteria: the qualitative description of the situations that the algorithm can handle (missing points, entry of points, etc), and the quantitative criteria given by the number of real structures found in the sequence (eg. the number of real trajectories, of real links between points, etc.), as well as the precision and recall of the algorithm for these different structures, defined by precision = # of correct structures found # of structures found ,

and recall = # of correct structures found # of actual structures .

The precision allows to measure the number of false positives (more precisely, 1-precision is the proportion of false positives amoung found structures), while the recall is linked to the number of false negatives (1recall represents the proportion of false negatives among actual structures). It is important to realize that the analysis of an algorithm must be done by considering simultaneously the precision and recall (or equivalent variables), since varying a parameter or a threshold of an algorithm generally does not improve both quantities but sets a different trade-off between the two, resulting in a better recall and worse precision or vice-versa.

In some experiments, the presence of noise points limits the interest of the number of real (whole) trajectories found as a significant criterion, although it is widely used in the litterature. Indeed, a well-placed noise point can sometime better fit the trajectory than its "real" counterpart, thus giving a realistic and usable trajectory as output, yet one that will not be counted as a real trajectory. We therefore chose to generally use the number of correct links as a significant structure for the precision and recall criteria. A link is simply two points that appear consecutively on a trajectory (possibly separated by a hole). Thus, if a noise point better fits a trajectory than its "real" counterpart, we will only "miss" two correct links (that include the real point), and "create" two false links (that include the noise point). However, when using the number of correct links, we do not account for trajectory over-segmentation, under-segmentation, or mixing. More precisely, if we split a trajectory in half, or if we join two distinct trajectories, we will barely notice it from the point of view of the number of correct links criterion, but we would have noticed it using the number of correct trajectories criterion. The same is true for "mixed" trajectories: if two trajectories cross at a point in time, we might start by following trajectory A, and then either choose to continue with trajectory A or to "hop" on trajectory B. In the latter case, the number of correct links criterion will barely be affected, but the number of correct trajectories criterion would. This particular problem of crossing trajectories appears however to be difficult to solve properly, and would certainly requires a priori knowledge. We believe that once the trajectories have been detected, even if they have been mixed, a simple post-processing task might be sufficient to split crossing trajectories in part at the crossing points, and reconstruct the real trajectories using an a priori (having the trajectories bounce if we are following billard balls, or having them cross if we are looking at fishes in an aquarium).

For the qualitative criteria, ROADS is able to account for missing and spurious points, as well as points leaving and entering the scene. The NFA algorithms come in two flavors, one that allows for missing points, and the other that does not. The latter is used for computational reasons (it is much faster) in some of the following experiments. Both NFA algorithms allow spurious points, as well as points leaving and entering the scene.

Trajectories without holes

Variable number of spurious points experiment

First, we investigate the influence of noise (spurious points) on trajectory detection. We generate sequences spanning 20 frames and having 20 trajectories, and we add 0 to 320 noise points uniformly at random to each image. Since we do not remove points, we can use the version of the NFA algorithm that does not take holes into account (Section 2). Then we run the NFA (ε = 1) and the ROADS (modes A,B,C) algorithms, and compare the results by computing the average recall and precision over 400 repetitions. For the precision estimate, the averaging is limited to the repetitions that lead to at least one detection, since the precision is not well defined when no structure is found.

As we explained earlier, the precision and recall are computed for two different criteria: the number of correct links and the number of correct trajectories. For the criterion based on the number of correct links (Fig. 10), the NFA algorithm performs much better in terms of precision: the precision remains very high (above 80%) for the NFA algorithm, but drops very fast for ROADS (under 20% when the number of spurious points exceeds 140). This illustrates a classical property of a-contrario detection models: the robustness to noise. As concerns the the NFA algorithm performs better than all versions of ROADS up to 200 spurious points, and is slightly under the C mode of ROADS beyond this level of noise. Considering the number of false detections made by ROADS at these levels of noise (the precision is under 10%), the global comparison remains in favor of the NFA algorithm. Table 4 clearly illustrates this: if ROADS manages to find a lot of correct links, it is solely because it makes a huge number of detections when the number of noise points increases, whereas the NFA algorithm correctly finds 20 trajectories in low noise and makes fewer detections when the noise level increases.

When we look at the number of correct trajectories found (Fig. 11), we see that ROADS is very good when there are no noise points. The NFA algorithm is a bit less efficient (both in terms of precision and recall) when the number of spurious points is under 40, but for higher levels of noise it is much more robust than ROADS, whose performances collapse very quickly (both in terms of precision and recall). A comparison with Fig. 10 is interesting here, because it shows that the "number of correct trajectories" is a very specific criterion that is not very relevant for complex or noisy data (the performances drop very quickly, much before the "number of correct links" is significantly affected). As we remarked before, there are plenty of reasons why a detected trajectory could be counted as undetected while it is very near an actual trajectory (a missing endpoint, a noise point fitting better the trajectory smoothness, trajectory crossings, NFA 20.1 20.2 18.9 15.5 7.1 6.1 ROADS(A) 20.0 70.0 151.5 227.1 300.9 335.7

Table 4 Average number of detected trajectories depending on the level of noise. We compare the average number of trajectories detected by the NFA and ROADS algorithms on data made of 20 real trajectories spanning the entire sequence plus a varying number of spurious points (from 0 to 320) in each frame. We see that NFA is very conservative in its detections (it only detects the trajectories that it considers to be non-ambiguous), and this results in a high precision (see Fig. 10). On the other hand, ROADS makes many false detections (it should not find more than 20 trajectories per sequence). etc.). Also, it is clear that applications based on data corrupted by a medium or high level of noise are more interested in a high rate of local point tracks (links) than in the perfect reconstruction of a very small proportion of the complete actual trajectories. This is why we shall not use the "number of correct trajectories" criterion any more in the following, but focus instead on the broader "number of correct links" criterion.

Variable density experiment

We now test how the algorithms behave when we increase the number of points. In this experiment, we do not consider spurious or missing points, so there is no noise and the difficulty of the trajectory detection problem only comes from the ambiguities produced by the large number of mixed trajectory points. We generate sequences of 20 frames, containing 10 to 140 points moving according to the PSMG model, where we allow trajectories to leave the image frame (when a trajectory leaves the image frame, we generate a new trajectory starting at a random position in the image frame, in order to keep a constant number of points throughout the sequence). Then we compute the precision and recall for the correct links criterion (Fig. 12) for the ROADS and NFA (without holes) algorithms.

When using the standard threshold (ε = 1) in the NFA algorithm, we obtain results which are similar (slightly better) than ROADS in terms of precision but significantly worse in terms of recall. However, knowing that there is no noise in these data, it makes sense to try to set the NFA threshold to +∞ (that is, no threshold), and in this case the results obtained by the NFA algorithms are similar to the best modes of ROADS. This is an unexpected good surprise for the NFA algorithm, which detects trajectories in a greedy way (by iterating a best-trajectory-detection/trajectory-removal process) without considering at all the global inter-frame assignment problem like ROADS. In the absence of noise points, one could have expected this assignment step to bring a significant edge to ROADS. On synthetic data containing 20 real trajectories spanning the entire sequence (20 frames) plus a varying number of spurious points (from 0 to 320), we compute the average recall (left) and precision (middle) obtained with the NFA and ROADS algorithms over 400 realizations, as a function of the level of noise (number of spurious points), or together (right). The most striking result here is that the precision of the NFA algorithm is almost constant, no matter the number of spurious points. This means that the NFA algorithm makes very few false detections (which is how we designed it), while keeping a recall rate that is above the one of ROADS as long as the number of spurious points is under 200 (which is more surprinsing). On the contrary, the poor precision of the ROADS algorithm in medium or high noise conditions makes its recall values quite unsignificant : if ROADS finds a large number of correct links, it is moslty because it proposes a high number of links, most of which are false detections (see also Table 4). 10 are now analyzed with a different criterion (the number of correct trajectories, instead of the number of correct links) for the definition of precision and recall. We can see that the NFA algorithm behaves much better than all ROADS modes as soon as there is a reasonable level of noise, while ROAD gives slightly better results when the noise level is very low. Note that the "number of correct trajectories" is a very specific criterion that is quite sensitive to local ambiguities (and to the selection of trajectory endpoints), and for that reason it is not very relevant when dealing with complex and/or noisy data. 12 Influence of the number of trajectories. The average recall (left) and precision (middle) are computed for the number of correct links criterion on 400 repetitions of synthetic data made of a given number of random trajectories (varying from 10 to 140) in a sequence of 20 frames. The analyzed algorithms are the three modes of ROADS (A, B, C), and two variants of the proposed NFA algorithm: the standard variant (threshold ε = 1 on the expected number of false alarms), and the no-threshold variant (ε = +∞). As we can see, the precision of both NFA variants is very high (like for ROADS B and C), but the recall of the standard NFA algorithm is significantly worse than the one of ROADS. In this setting where no noise points are present, these missing detections can be avoided by removing the thresholding process in the NFA algorithm: for this ε = +∞ variant, both recall and precision are as good as the best modes of ROADS.

Sensitivity to data smoothness

The trajectories generated in the previous experiments are somehow smooth and quasi-linear (see Fig. 9). In order to see if the algorithms can cope with trajectories that do not completely fit the model, we try to detect trajectories having a potentially high acceleration. For that purpose, we consider different values of σ , the standard deviation of the acceleration magnitude used in the PSMG synthesis procedure (see the very begining of Section 4.2). The effect of this parameter on the synthetized trajectories is illustrated in Fig. 13.

We first reproduce the last experiment (Fig. 12), in which there are no noise points and the number of synthetized trajectories ranges from 10 to 140, and evaluate the effect of the σ parameter both for the unthresholded NFA algorithm (Fig. 14) and the ROADS B algorithm (Fig. 15). Whereas the performance of the NFA algorithm barely depends on σ (and remains high), ROADS exihibits a high sensitivity to this parameter, and its performance quickly collapses as σ increases. The same conclusion arises from the analysis of data generated with 10 noise points per frame plus 20 synthetic trajectories (see Fig. 16).

Thus, the sensitivity to data smoothness is a major difference between the NFA and ROADS algorithms. The poor results obtained by ROADS for σ = 1 (see Fig. 15) could probably be improved by a specific choice of the ROADS parameters (specially adapted to σ = 1), but this kind of optimization will not be efficient on most real-world data, since one can expect to observe a high variability of accelerations on such data. Conversely, the robustness of the NFA algorithm to the σ parameter is an indication that it can probably handle well real-world data containing various levels of acceleration.

Parameter tuning

One major interest of most a-contrario methods is that they yield "parameterless" detection algorithms, or, more pre-Fig. 13 Changing the acceleration variance. A sample of 20 trajectories generated using the PSMG algorithm, when the standard deviation of the acceleration magnitude (σ ) is 1 (left) and 4 (right). The points of two trajectories have been highlighted to show the speed. We study the sensitivity of the algorithms to data variability by analyzing their performances when we increase σ . cisely, algorithms for which there exist natural values of the parameters that work well in all situations. Both NFA algorithms we propose here (the no-hole and hole versions) have only one parameter: the threshold ε used to decide whether a trajectory should be detected or not. Since ε corresponds to an upper bound on the expected number of false alarms in pure noise data, its default value is classically set to 1 (see [START_REF] Desolneux | From Gestalt Theory to Image Analysis: A Probabilistic Approach[END_REF]). In Fig. 17, we examine the sensitivity of the NFA no-hole algorithm with respect to the choice of ε. We use the same experimental setting as in Fig. 10 (20 frames containing 20 real trajectories plus several spurious points), and examine how recall and precision are affected by different choices of ε. The results clearly show that the default value ε = 1 (log 10 ε = 0) is nearly optimal, in the sense that it is small enough to guarantee a strong precision control, and large enough to offer good recall performances. It is nonetheless interesting to notice that slightly better performances (same precision and better recall) can be obtained with greater values of ε (typically log 10 ε = 2 or 3).

In Fig. 18, the average precision/recall curve obtained with the NFA algorithm for different values of the threshold log 10 ε (in the case of 160 spurious points) is displayed on the left. On the right, we report the average performances of ROADS on the same data points, considered for several values of the two main parameters of this algorithm, namely the maximal speed and the maximal smoothness. The maximum speed parameter varies from the actual value in the [-50%, +50%] range (from one curve to another) and the maximum smoothness varies from the actual value in the [-95%, 50%] range (inside each curve). Note that the best performances of ROADS (-25% speed, -90% smoothness) are obtained inside these ranges.

As we can see, not only the performances of ROADS are way under those of NFA on these data, but also the parameter tuning is much more difficult and crucial (we have to explore carefully a bidimensionnal domain, while ε = 1 is almost optimal for the NFA algorithm).

NFA as a criterion for trajection selection

Contrary to ROADS, which is by nature an algorithm (relying in particular on some heuritics), the NFA we propose here is first and foremost a criterion to compare trajectories. The greedy algorithm we described, based on the iteration of a "best trajectory (minimal NFA) detection / trajectory removal" process, is only one possibility to use the NFA criteria ( 5) and ( 23), and it is possible to design other algorithms based on these criteria. In particular, given a trajectory detection algorithm, it is always possible to use the NFA criterion as a post-processing step that simply keeps from the output of the considered algorithm the trajectories having a NFA under a certain threshold ε. 19 ROADS output filtered by the NFA algorithm. We consider the same synthetic data as in Fig. 18, but now add to the comparison of NFA and ROADS algorithms a combination of them that consists in detecting trajectories with ROADS and keeping only those having a NFA under a certain threshold ε. Since each algorithm depends on parameters (1 for NFA, 2 for ROADS, 3 for ROADS+NFA), we explore systematically all parameter values and compute the upper performance enveloppe (curves named best). As we can observe, the major drawback of ROADS (which is its high rate of false detections) can be corrected by NFA fitering, which results in a dramatic increase of precision (up to the level of the NFA algorithm alone). However, this correction does not permit to attain the same level of recall (around 0.75 for NFA, versus 0.6 for ROADS+NFA in the high precision zone). Note also that the mixed ROADS+NFA algorithm would be much more complicated to use than NFA alone, due to the 3 parameters that have to be set. We generate 20 trajectories spanning the whole sequence (20 frames), and remove randomly 20% percent of the points, before we add a varying number of spurious points (from 0 to 70). On these synthetic data (with 400 repetitions), we estimate the recall (left) and the precision (middle) of the ROADS and NFA algorithms for the "number of correct links" criterion. The obtained results are very similar to those of Fig. 10: the recall values are roughly the same for all algorithms, but only the NFA algorithm manages to maintain a high precision (above 0.9) as the number of spurious points increases, while all ROADS variants make lots of false detections.

Trajectories of real-world images

The snow sequence

In this part, we evaluate the relative performances of NFA and ROADS algorithms on a real-world sequence named snow. To produce this sequence, we filmed falling snowflakes in front of a dark metal door with a high-speed (210 frames per second) camera, and then subsampled the high speed sequence at 30 fps by taking 1/7 of the original frames. This way, we obtained a classical 30 fps sequence of 40 images of 480 × 360 pixels, on which we ran a simple point extraction process that we describe below. The high-speed version was processed in the same way and used in order to build a hand-made ground truth for trajectories.

We purposefully used a very simple extraction process to produce data as objectively as possible, without trying to adapt the detection algorithm in a way that would affect (and ease) the tracking part. The snowflakes (but also some stains on the metal door background) were detected in the following way: we smoothed the images using a simple Gaussian kernel, and we computed the mean background image on a few frames of the subsampled (30 fps) sequence. We then thresholded the image differences, processed the result with a morphological closing, and extracted the connected components. For each connected component, we kept the centroid position, rounded to the nearest point on the integer grid, as a trajectory data point.

In the resulting point sequence, many objects were detected as several close points in the sequence (in particular the stains on the background and some big snowflakes). This made it sometimes hard to define the ground truth trajectories. To alleviate this difficulty, we kept only isolated points by removing all points in the sequence that were below a certain distance of another point (we chose the smallest radius that would resolve almost all ambiguities, in practice 12 pixels). An example of detections on one frame of the sequence is displayed in Fig. 21. We finally extracted the ground truth trajectories by hand.

The resulting point sequence is interesting because it presents a mix of difficulties: there are widely varying trajectory types (points in the background which practically do not move, very slow snowflakes with curvy trajectories, very fast snowflakes with almost linear trajectories). There are missing points (missing detections or detections removed because of the non-isolated point removal process), and a few noise points (but the relatively high detection threshold gave more missing points and fewer noise points).

Finally, we subsampled the high-speed point sequence by keeping only 1/7 of the frames, and subsampled accordingly the ground truth trajectories. The resulting trajectories containing less than 3 points were eliminated from the ground truth reference, but the corresponding points were kept in the data (thus becoming noise points). The final result of this process (30 fps snow point sequence and associated ground truth) is available on the website http://www.mi.parisdescartes.fr/ ~moisan/astre/ The first row of Fig. 25 gives an idea of the ground truth trajectories extracted from the snow sequence.

Parameter tuning

There are several parameters to set for ROADS (see Table 2), and different settings lead to varying results. Namely, we can set the size s of the time scope (we chose 2, giving the best results), the minimal number p min of consecutive present points for a trajectory to be considered (we chose 1, 3, 5 or 7), the maximal length of interpolation a max before we loose the trajectory (we chose 0, 4, 8 or +∞), the maximal smoothness criterion φ max and the maximal speed d max . The way to choose the best parameters is not obvious, but it appears in Fig. 22 that the most important parameter is the maximum allowed speed d max . The choices p min = 1, a max = 0 and φ max = 0.4 are among the best possible for the snow point sequence, and would probably achieve reasonable performances on similar sequences too. As concerns the choice of d max , the ground truth value ( 160) is much too large, and much better results are obtained with d max = 15. This fact, which comes from the unability of ROADS to deal with a variety of trajectory speeds at the same time, is analyzed more precisely later. Note that the ground truth value of φ max is 0.58.

On the snow sequence, extracting trajectories using the two NFA algorithms (the one with holes and the one without holes) would return trajectories having a value of log 10 (NFA) varying from -40 to +10, and the "optimal" precision/recall values would be obtained (for both algorithms) by thresholding this value with log 10 ε = +5 (see Fig. 23). Even without access to the ground truth, finding this value is relatively easy, since one simply has to look for values slightly above the (nearly optimal) default value log 10 ε = 0. This strategy works well in all synthetic experiments we considered earlier, and also in the present case of the snow sequence. In view of the false detection control offered when log 10 ε = 0, such a strategy is probably efficient on most (if not all) point sequences.

Thus, as we mentioned before, one great interest of the NFA algorithm is that the parameter tuning step is much more easier than in other algorithms like ROADS, for which it can be a real burden, especially when dealing with complex data (with unknown ground truth) on which the effect of a parameter change can be very difficult to evaluate. This relative parameter sensitivity is illustrated in Fig. 24.

Comparison of ROADS and NFA algorithms

To compare the results obtained by the ROADS and NFA algorithms on the snow sequence, we use for each algorithm two different settings: the default setting and the best setting.

For ROADS, the default setting corresponds to a max = +∞, p min = 1, d max = 130, and φ max = 0.58. Note that φ max = 0.58 corresponds to the oracle value, that is, the (theoretically unknown) maximum value of φ on the ground truth trajectories. For d max , we chose the value d max = 130 to allow ROADS to detect all the trajectories in the main bulk of trajectories (choosing d max as the real maximal speed (160) would give worse results). The best setting for ROADS was chosen after a careful (and a bit cumbersome) parameter analysis (see Fig. 22), which leads to a max = 0, p min = 1, φ max = 0.6, and d max = 20.

Concerning the NFA algorithms, the default and best settings simply correspond to log 10 ε = 0 and log 10 ε = +5 re-spectively (see Fig. 23). As for ROADS, the best setting was chosen a bit arbitrarily, as one of the best recall values avoiding a significant loss of precision. The strong L-shaped aspect of the precision-recall curves made that choice quite easy (in the sense that other possible choices would not differ much).

The results obtained with ROADS and NFA algorithms on the snow sequence (for both default and best settings) are shown in Fig. 25, and the precision-recall performances reported in Table 5. In the best configuration, all algorithms attain similar performances in terms of precision, but at the price of a high number of false detections (poor recall) for the ROADS algorithm. Moreover, the default parameter setting of ROADS gives very poor results, while the default parameter setting of NFA (ε = 1) achieves a better performance than the best ROAD settings. These results are analyzed in greater details in Fig. 26, where it appears that the main limiting factor of the ROADS algorithm seems to be its unability to handle simultaneously (that is, with a same set of parameters) trajectories with various lengths and speeds. The ROADS and NFA algorithms are run on the snow sequence with their default and best settings, and their performances are analyzed in terms of recall and precision. In order to permit a more accurate analysis, separate recall scores are also computed by considering separately fast foreground (fg) objects (snowflakes) and the almost static background (bg) objects (staints on the background door). As we can see, the comparison is clearly in favor of the NFA algorithm. With the best settings, the obtained precision is roughly the same, but the ROADS algorithm is unable to achieve an interesting detection rate on the foreground objects (the snowflakes), which results in a poor overall recall.

We end this section with a little discussion on computation times. On the snow sequence, the ROADS algorithms runs in about 10 seconds for the default parameter values, and in 0.1 second for the best parameter values (this does not take into account, of course, the time needed to find these best parameters values). In comparison, the NFA "no-hole" algorithm runs in 1.1 second on the same data, while the NFA "hole" algorithm takes 35 minutes. The speed ratio is large, but it should be noted that in several applications such a computation time is not a problem, because the production of the point sequence may take much more time (consider a biological experiment relying on cell tracking for example). Moreover, there is an intermediate way of speeding up the algorithm while keeping the interesting property of allowing holes (which may be crucial for some data). It consists in using the NFA with holes criterion (Equation 23) while limiting the exploration of trajectories to those that do not have holes longer than h (this requires only a very simple modification of the algorithm). In this modified algorithm, the integer parameter h has no influence on the detection thresholds, but greatly decreases the computation time by limiting the search to the most common trajectories. It may be very useful to find a tradeoff in situations where computation time matters more than full detection performances, or in a quick-analysis stage used before running the exact (h = ∞) NFA "hole" algorithm. On Table 6, we reported the influence of the h parameter on the computation time and detection performances (precision/recall) for the snow sequence (log 10 ε = 5). 6 Limiting the size of holes in the NFA algorithm. This table reports the computation time and the precision/recall performances obtained after running on the snow sequence the NFA "hole" algorithm (log 10 ε = 5) restricted to trajectories having no hole longer than h (for h = 0, the "no-hole" NFA algorithm is used). As h decreases, the computation time drops quickly, while maintaining a high level of precision and recall on that sequence.

Conclusion

We presented two point-tracking algorithms based on the acontrario framework, which are able to detect trajectories in point sequences without additional information. The first algorithm (NFA "no-hole") is restricted to complete trajectories (that is, without holes), while the second (NFA "hole") can recover trajectories with missing points. Both algorithm are very robust to noise, in the sense that they are designed to avoid halucinating trajectories in noise data. Another strength of these algorithms, which comes from the a-contrario approach, is that they do not require to set parameters: even the only threshold that can be set (ε, which balances between precision and recall) may be left to its default value (ε = 1). When compared to the state-of-the-art ROADS algorithm, these two algorithms perform very well, both on simulated and real-world data. In particular, they show a very high level of precision (that is, a very low rate of false alarms) while maintaining a good level of recall (actual trajectories detected). The absence of required parameter setting and the robustness to trajectory variability (speed, length, acceler- First row: trajectories from ground truth (right) and a sample of ground truth trajectories where successive points are represented as dots to give an idea of the objects speeds. Second row: trajectories found by the NFA algorithm, for the default parameter value (ε = 1, left) and for the optimal parameter value (log 10 ε = 5, right). The well-detected trajectory links are drawn in black, and the wrong ones in gray. Third row: trajectories found by the ROADS algorithm, for the default parameter values (left) and the optimal parameter values (right). As we can observe, the ROADS algorithm makes many false detections with the default parameter values, and very few detections with the optimal parameter values. On the contrary, the NFA algorithm with the default parameter value finds a large part of the actual trajectories and yields almost no false detection (see Fig. 23). When using the best parameter value, almost all trajectory links are found, and only a few spurious detections occur.

ation) appears to be a very useful feature when real-world point sequences have to be processed.

Compared to several classical tracking algorithm, the two NFA algorithms we proposed have the advantage that they are based on the exact minimization (and thresholding) of a simple criterion, without any heuristic or approximation. This criterion not only permits to rate the quality of each trajectory, but might also be used in a very simple way to remove false detections from to the output of another algorithm.

The "no-hole" algorithm is rather fast on standard data, but the principal limitation of the "hole" algorithm is its high computational and memory costs. However, these costs may be mitigated by using the easily parallelizable nature of the 5) are analyzed in function of the maximal trajectory speed, for both best (left) and default parameter settings.

The possible values of the maximal trajectory speed are divided into bins (horizontal axis), and the blue histogramms indicate the number of corresponding actual trajectories. Then, the recall of each algorithm is analyzed inside each bin, in red for NFA and in dashed green for ROADS.

As we can see, the ROADS algorithm does not manage to handle simultaneously trajectories with various speeds: the detection is focused either on trajectories with middle-range speeds (default setting), or on very slow trajectories (best setting). The NFA algorithm, which combines the trajectory smoothness and length into a single NFA criterion (hence avoiding a speed threshold), does not suffer from this dilemma, as it clearly appears on the left (best setting) graph.

algorithm, or by reducing the size of the trajectory search space with additional constraints (for example, a bound on the maximum number of consecutive missing points). We chose in this paper to use a trajectory smoothness criterion based on the maximum acceleration, but similar approaches could be developped as well for other local criteria (speed, angular acceleration, etc.) and for other global costs (a sum cost instead of a max cost for example). The presentation was made in a two-dimensional setting, but the generalization to higher dimensions (3D points, or more) is straightforward. Note also that the proposed framework could probably be extended to the case when points come with features (intensity, shape, etc.), in the same spirit as [START_REF] Noury | Determining point correspondences between two views under geometric constraint and photometric consistency[END_REF] extended the framework of [START_REF] Moisan | A probabilitic criterion to detect rigid point matches between two images and estimate the fundamental matrix[END_REF] for Fundamental Matrix Estimation.

Fig. 3

 3 Fig. 3 Asymptotic and non-asymptotic critical accelerations. The critical acceleration δ c (here, in base-10 log scale) is an increasing function of ℓ that approaches its upper bound (dotted line) ∆ = 1 √ Nπ when ℓ = K. Here K = 50 and the three curves correspond to N = 15, N = 50, and N = 200 respectively.

Fig. 5

 5 Fig.5Minimal number of frames required for detectability. Theses three curves represent, as a function of the number of points N, the minimum number of frames (K) required for a trajectory with length ℓ = K and various values of the maximal acceleration δ to be 1-meaningful according to the NFA criterion (15).

Fig. 6

 6 Fig. 6 Comparison between the NFA and NFA ′ models(15, 16). Each curve represents, as a function of µ, the maximal value of α = Nπδ 2 allowed for a trajectory with length ℓ = ⌊µK⌋ and maximal acceleration δ to be 1-meaningful in a sequence of K images, among N = 100 points per frame. The red upper curves are obtained with the criterion NFA, whereas the blue lower curves are obtained with the criterion NFA ′ . The full curves correspond to K = 100, and the dashed curves correspond to the asymptotical estimates obtained when K → +∞, that is, log α = 0 for NFA and log α = log N • (µ -1)/µ for NFA ′ . These curves clearly demonstrate that not only NFA is better suited to the detection of small trajectories than NFA ′ (it allows for trajectories having a much larger maximal acceleration), but it is also more efficient even for relatively large values of µ (NFA ′ being slightly better only for almost complete trajectories). Asymptotically, NFA is always the best choice.

  Fig.10Influence of spurious points (# correct links criterion). On synthetic data containing 20 real trajectories spanning the entire sequence (20 frames) plus a varying number of spurious points (from 0 to 320), we compute the average recall (left) and precision (middle) obtained with the NFA and ROADS algorithms over 400 realizations, as a function of the level of noise (number of spurious points), or together (right). The most striking result here is that the precision of the NFA algorithm is almost constant, no matter the number of spurious points. This means that the NFA algorithm makes very few false detections (which is how we designed it), while keeping a recall rate that is above the one of ROADS as long as the number of spurious points is under 200 (which is more surprinsing). On the contrary, the poor precision of the ROADS algorithm in medium or high noise conditions makes its recall values quite unsignificant : if ROADS finds a large number of correct links, it is moslty because it proposes a high number of links, most of which are false detections (see also Table4).

  Fig.11Influence of spurious points (# correct trajectories criterion). The results of the experiments conducted in Fig.10are now analyzed with a different criterion (the number of correct trajectories, instead of the number of correct links) for the definition of precision and recall. We can see that the NFA algorithm behaves much better than all ROADS modes as soon as there is a reasonable level of noise, while ROAD gives slightly better results when the noise level is very low. Note that the "number of correct trajectories" is a very specific criterion that is quite sensitive to local ambiguities (and to the selection of trajectory endpoints), and for that reason it is not very relevant when dealing with complex and/or noisy data.

  Fig.12Influence of the number of trajectories. The average recall (left) and precision (middle) are computed for the number of correct links criterion on 400 repetitions of synthetic data made of a given number of random trajectories (varying from 10 to 140) in a sequence of 20 frames. The analyzed algorithms are the three modes of ROADS (A, B, C), and two variants of the proposed NFA algorithm: the standard variant (threshold ε = 1 on the expected number of false alarms), and the no-threshold variant (ε = +∞). As we can see, the precision of both NFA variants is very high (like for ROADS B and C), but the recall of the standard NFA algorithm is significantly worse than the one of ROADS. In this setting where no noise points are present, these missing detections can be avoided by removing the thresholding process in the NFA algorithm: for this ε = +∞ variant, both recall and precision are as good as the best modes of ROADS.

  Fig.19ROADS output filtered by the NFA algorithm. We consider the same synthetic data as in Fig.18, but now add to the comparison of NFA and ROADS algorithms a combination of them that consists in detecting trajectories with ROADS and keeping only those having a NFA under a certain threshold ε. Since each algorithm depends on parameters (1 for NFA, 2 for ROADS, 3 for ROADS+NFA), we explore systematically all parameter values and compute the upper performance enveloppe (curves named best). As we can observe, the major drawback of ROADS (which is its high rate of false detections) can be corrected by NFA fitering, which results in a dramatic increase of precision (up to the level of the NFA algorithm alone). However, this correction does not permit to attain the same level of recall (around 0.75 for NFA, versus 0.6 for ROADS+NFA in the high precision zone). Note also that the mixed ROADS+NFA algorithm would be much more complicated to use than NFA alone, due to the 3 parameters that have to be set.

  Fig.20Influence of spurious points for trajectories with holes. We generate 20 trajectories spanning the whole sequence (20 frames), and remove randomly 20% percent of the points, before we add a varying number of spurious points (from 0 to 70). On these synthetic data (with 400 repetitions), we estimate the recall (left) and the precision (middle) of the ROADS and NFA algorithms for the "number of correct links" criterion. The obtained results are very similar to those of Fig.10: the recall values are roughly the same for all algorithms, but only the NFA algorithm manages to maintain a high precision (above 0.9) as the number of spurious points increases, while all ROADS variants make lots of false detections.

Fig. 21

 21 Fig. 21 An image of the snow sequence (inverted grayscale), with overlayed detections.

  Fig.25Trajectories found by ROADS and NFA algorithms on the snow sequence. First row: trajectories from ground truth (right) and a sample of ground truth trajectories where successive points are represented as dots to give an idea of the objects speeds. Second row: trajectories found by the NFA algorithm, for the default parameter value (ε = 1, left) and for the optimal parameter value (log 10 ε = 5, right). The well-detected trajectory links are drawn in black, and the wrong ones in gray. Third row: trajectories found by the ROADS algorithm, for the default parameter values (left) and the optimal parameter values (right). As we can observe, the ROADS algorithm makes many false detections with the default parameter values, and very few detections with the optimal parameter values. On the contrary, the NFA algorithm with the default parameter value finds a large part of the actual trajectories and yields almost no false detection (see Fig.23). When using the best parameter value, almost all trajectory links are found, and only a few spurious detections occur.

  Fig.26Detailed analysis of ROAD and NFA recall performances. The recall performances of ROADS and NFA "hole" algorithms on the snow sequence (as given in Table5) are analyzed in function of the maximal trajectory speed, for both best (left) and default parameter settings. The possible values of the maximal trajectory speed are divided into bins (horizontal axis), and the blue histogramms indicate the number of corresponding actual trajectories. Then, the recall of each algorithm is analyzed inside each bin, in red for NFA and in dashed green for ROADS. As we can see, the ROADS algorithm does not manage to handle simultaneously trajectories with various speeds: the detection is focused either on trajectories with middle-range speeds (default setting), or on very slow trajectories (best setting). The NFA algorithm, which combines the trajectory smoothness and length into a single NFA criterion (hence avoiding a speed threshold), does not suffer from this dilemma, as it clearly appears on the left (best setting) graph.

Table 2

 2 Parameters used in the ROADS algorithm.

		smoothness model parameter
	d max	maximal allowed speed
	φ max	maximal allowed smoothness criterion
	s	scope width parameter
	a max	max. # of missing consecutive points on a track
	p min F γ g , F l γ	min. # of present consecutive points on a track optimization cut-off constants
	h max	max. # of hypotheses made when optimizing

Table 5

 5 Performances of ROADS and NFA algorithms on the snow sequence.

	mode	algorithm	recall (fg, bg)	precision
	default	ROADS NFA no holes 0.50 (0.44, 0.55) 0.09 (0.19, 0.00)	0.08 0.99
		NFA holes	0.55 (0.44, 0.64)	0.91
	best	ROADS NFA no holes 0.71 (0.74, 0.68) 0.40 (0.13, 0.65)	0.96 0.95
		NFA holes	0.76 (0.78, 0.74)	0.92

Acknowledgments

We are very grateful to C.J. Veenman for providing us with a C implementation of the ROADS algorithm.

This work was funded by the French National Research Agency (ANR) under contract ANR-09-PIRI-0030-03 (PAGDEG project).

Fig. 14 Influence of the data smoothness (unthresholded NFA algorithm). We use the same experimental setting as Fig. 12, and examine the sensitivity of the unthresholded (ε = ∞) NFA algorithm to the smoothness of the analyzed synthetic data. More precisely, we consider several values of σ , the standard deviation of the acceleration magnitude (a parameter of the synthesis algorithm, PSMG), and estimate the precision and recall (correct links criterion) as functions of the number of synthetic trajectories. We can see that the NFA algorithm is extremly robust to σ , since both the precision and recall performance curves remain unchanged when σ varies. We analyze the same data as in Fig. 14, now with the ROADS algorithm, mode B (the best mode for these data, see Fig. 12). Contrary to what happens for the NFA algorithm, the ROADS method exhibits a severe sensitivity to σ , since both recall and precision performances, which were at the same level as the NFA algorithm for σ = 0.2 (grey shadow curves), are strongly affected when σ increases. As we shall see in Section 4.5, the sensitivity/robustness to data variability has strong consequences when real-world data are analyzed. Note incidentally the strong similarity between the recall and the precision curves, which comes from the fact that in this set of experiments, the number of detected links is most of the time equal to the number of actual links (see Equation ( 31) and (32)), probably because there are no spurious points. 16 Sensitiviy to data smoothness. The average recall and precision of the ROADS algorithm (modes A, B, C) and the standard NFA algorithm (ε = 1) are compared in synthetic data made of 10 noise points per frame plus 20 random trajectories spanning 20 frames, in function of the standard deviation of the acceleration magnitude (σ ). These results corroborate the ones obtained in Fig. 14 and 15: the performances of the NFA algorithm are not too much affected by the increase of σ (except for the recall when the variance becomes large, probably because the problem of recovering the true trajectories becomes objectively difficult), whereas the performances of all ROADS algorithms collapse, both in terms of precision and recall. Fig. 17 Influence of the NFA threshold (ε). We consider the same experiment as in Fig. 10 (that is, 20 real trajectories spanning 20 frames, with a given number of spurious points in each frame), and examine the influence of the threshold ε arising in the NFA algorithm. Recall and precision curves are plotted in function of the number of spurious points, for different values of log 10 ε (ranging from -4 to +∞). We can see that the good precision control predicted by the theory for ε 1 (log 10 ε 0) is well achieved, since the first significant precision losses occur around log 10 ε = 3. Hence, the default value log 10 ε = 0 is a good compromise in this experiment, even if slightly better recalls (without significant precision losses) can be achieved by using greater values like log 10 ε = 2. 18 Performance and parameter tuning. We consider a particular case of Fig. 17, that is, synthetic data made of sequences of 20 frames containing 20 real sequences and 160 spurious points on each frame. The average performances in terms of precision/recall is then evaluated for the NFA algorithm (left) and the ROADS algorithm (right), with varying values of the algorithm parameters. For the NFA algorithm, the only parameter is the treshold ε (or log 10 ε, as displayed on the figure), and we can see that the default value log 10 ε = 0 is very near to be optimal, as was remarked earlier in the comment of Fig. 17. For ROADS, not only the performances are lower (especially in terms of precision), but they are also quite sensitive to the choice of the maximum speed and maximal smoothness parameters.

We tested this possibility with the ROADS algorithm, and reported in Fig. 19 the results obtained on the synthetic data used in the previous section (parameter tuning). It appears that the mixed ROADS+NFA algorithm we obtain this way performs much better than ROADS alone in terms of precision (because the NFA filtering permits to eliminate most false detections), but the performances in terms of recall do not attain the ones of the NFA algorithm alone. Hence, the "NFA filtering" strategy is efficient but does not provide a particularly interesting new algorithm when applied to ROADS. It is not impossible, however, that such a strategy could be successful, in particular in situations where only special kinds of trajectories appear and a good detection algorithm (in terms of recall) exists. In that kind of situation, one could expect NFA filtering to increase the precision up to a high level, without damaging to much the recall performances. Note that such a strategy guarantees, thanks to the properties of the NFA criterion (1), the control of the number of false detections in random data.

Trajectories with holes

We now examine the performances of the second NFA algorithm (Section 3), which is able to handle trajectories with holes. We compare it to ROADS using the same kind of conditions as in Fig. 10 (20 real trajectories, 20 frame, several spurious points added in each frame), except that we now consider uncomplete trajectories (20% of the points of the true trajectories are removed before spurious points are added). The conclusions made in the no-hole case remain unchanged (see Fig. 20) : the ROADS algorithm detects true trajectory links as well as the NFA algorithm, but at the price of many false detections, whereas the NFA algorithm makes almost no false detection (the precision remains above 0.9, even for 70 spurious points per frame). We vary all ROADS parameters on the snow point sequence, and show the associated performances in the (recall,precision) plane using the available ground truth for that sequence. Each column has a distinct a max = 0, 4, 8, +∞, and each row has a distinct p min = 1, 3, 5, 7. Each curve corresponds to a different maximal smoothness criterion value φ max = 0.2, 0.4, 0.6 and each point of a given curve corresponds to a different maximal speed criterion d max = 2, 5, 10, 15, 20, 25, 30, 40, 80. The big red dot corresponds to the parameters φ max = 0.4 and d max = 15, which seem to achieve a good precision/recall compromise for all values of a max and p min . The numbers indicate the corresponding recall and precision. 23 NFA parameter tuning on the snow sequence. The performances of the two NFA algorithms (with holes and without holes) on the snow sequence are represented in the (recall,precision) plane in function of the threshold parameter log 10 ε. While the precision remains merely constant, a good recall is obtained with the default value (0) of log 10 ε, but the results can be improved by choosing a slightly greater value (log 10 ε = +5), which corresponds for the two algorithms to the "NFA best" point. 24 Comparison of NFA and ROADS algorithms on the snow sequence. All results obtained on the snow sequence with the NFA (with holes) and ROADS algorithms are represented in the (recall,precision) plane, with a point for each set of parameters (thus, the optimal performance of each algorithm is the curve obtained as the upper-right enveloppe of its points). We can see not only that ROADS is much more sensitive to the parameter choice than the NFA algorithms, but also that its overall performance in terms of recall is significantly worse, even with an optimal choice of its parameters.

Available resources

The snow point sequence, as well as open-source codes for the proposed algorithms are freely available on the ASTRE (A-contrario Smooth TRajectory Extraction) website http://www.mi.parisdescartes.fr/ ~moisan/astre/ It contains a simple Python implementation, an optimized one in C language working under Linux Operating System, and a detailed installation and usage manual.