Daniel Chen

Leonidas Guibas

Xiaoye Jiang

Christian Sommer

Data-Driven Trajectory Smoothing

Keywords: Algorithms, Certication, and Cryptography Algorithmics, Programming, Software and Architecture Geographic Information System (GIS), Theory, Smoothing, Trajectories, Data-Driven Techniques Geographic Information System (GIS), Theory, Smoothing, Trajectories, Data-Driven Techniques

Motivated by the increasing availability of large collections of noisy GPS traces, we present a new data-driven framework for smoothing trajectory data. The framework, which can be viewed of as a generalization of the classical moving average technique, naturally leads to ecient algorithms for various smoothing objectives. We analyze an algorithm based on this framework and provide connections to previous smoothing techniques. We implement a variation of the algorithm to smooth an entire collection of trajectories and show that it perform well on both synthetic data and massive collections of GPS traces.

Figure 1: The underlying path in solid black, a collection of traces in dashed red, and the query trajectory in dashed-dotted blue. Note that the noisy traces only partially sample the underlying path and that there may not exist a single trace sampling the same portion of the path as the query trajectory. 1 Introduction

Motivation Large scale collections of trajectory data, such as GPS traces, are becoming widely available, and many applications [START_REF] Zhang | Integration of gps traces with road map[END_REF][START_REF] Zheng | Learning travel recommendations from user-generated gps traces[END_REF][START_REF] Huang | Activity identication from gps trajectories using spatial temporal pois' attractiveness[END_REF][START_REF] Jeung | Discovery of convoys in trajectory databases[END_REF][START_REF] Chen | Discovering popular routes from trajectories[END_REF][START_REF] Yuan | T-drive: driving directions based on taxi trajectories[END_REF][START_REF] Yoon | Smart itinerary recommendation based on user-generated gps trajectories[END_REF] that involve understanding and extracting information out of such data are emerging.

However, using such data in practice is often challenging as the noise levels often exceed those for which many algorithms are designed. One solution to cope with this problem is a class of tools known as smoothing techniques, which preprocess the data by removing the noise, and thereby revealing the important underlying paths. Indeed, in practice, standard smoothing techniques, such as kernel smoothers [START_REF] Hastie | The Elements of Statistical Learning[END_REF], smoothing splines [START_REF] Hastie | Generalized additive models[END_REF], generalized additive models [START_REF] Hastie | Generalized additive models[END_REF],

and Kalman lters [START_REF] Welch | An introduction to the kalman lter[END_REF] are often used to preprocess GPS data in preparation for road network generation [START_REF] Bruntrup | Incremental map generation with GPS traces[END_REF] and other model extraction tasks [START_REF] Huang | Meta: A mobility model of metropolitan taxis extracted from gps traces[END_REF].

However, these methods do not leverage the special structure of trajectory data and the opportunities for statistical estimation given by the large collection of trajectories at the same time. In this paper, we explore a new smoothing framework specically tailored for trajectory data. Moreover, our approach is data-driven rather than model-driven, and hence is suitable for the setting where we have a large collection of traces sampling underlying paths which may belong to vastly dierent classes of curves. In this paper, we formalize the data-driven smoothing paradigm by developing both a model for the input data, and a novel and elegant algorithmic framework that results in practical and ecient algorithms.

Input and Output We assume that there exist: (1) An underlying path which we cannot observe directly, (2) a database of traces partially sampling the path and (3) a query trajectory sampling from the path that we seek to smooth using the traces (See Figure 1). The goal is then to use the collection of traces to output a smoothed version of the query trajectory. Note that there may not be a single trace that samples the entire portion of the path corresponding to the query trajectory. In the particular example of GPS data, this model allows the query trajectory to correspond to a route in the road network that we do not have in the database. However, the route can be split into many sub-routes, where each sub-route is covered by a trace in the database. If the particular application requires that the entire database of traces be smoothed, we can iterate over the traces, and assign each as the query trajectory.

RR n°7754

Algorithmic Framework Inspired by [START_REF] Takens | Detecting strange attractors in turbulence[END_REF], we propose a simple and elegant method that begins by embedding each point p of the traces and query trajectory into a high dimensional space as p using its delay coordinates, which are the coordinates of a window of sample points preceding and following p (see Section 3). We call this process lifting p, the high dimensional space the lifted space, and the embedded image p the lifted image of p. We then move each point p in the lifted image of the query trajectory towards several of its nearest neighbors in the lifted space (Laplacian smoothing). The exact set of nearest neighbors selected for this step depends on the smoothing objective. Finally, we recover the original traces by averaging the relevant coordinates of the lifted images.

Analysis and Experiments We prove that, under natural smoothing objectives, the process of using delay coordinates achieves much better results than considering the individual sample points by themselves. We show that under some weak probabilistic assumptions, a variant of the algorithm moves a noisy query trajectory towards the portion of the underlying path that it is sampled from. Under certain assumptions, the technique can be shown to be an unbiased estimator of the underlying path with high probability. Furthermore, we are able to dramatically reduce the variance of the error from that of the original data by smoothing with respect to a large database of traces.

We also show that our method can be viewed of as a generalization of standard moving average techniques for smoothing time series data. When there is no database of traces available to smooth a query trajectory, then our method is equivalent to the moving average technique. The result of such a smoothing technique optimizes an objective function that weighs the distance between the query trajectory and a constant velocity trajectory in a natural manner.

Experimentally, we consider both synthetic noisy data, as well as a collection of real GPS traces from Moscow and a collection of taxi cab traces from an unidentied major city. We adapt our framework for all three datasets and show that they result in trajectories that are signicantly smoother and less noisy than the input trajectories.

Related Work Our work is most similar to, and motivated by that of Cao and Krumm [START_REF] Cao | From gps traces to a routable road map[END_REF], who also smooth GPS traces before using them to construct a routable road map. Their approach, which is based on simulating forces of physical attraction, is also data-driven rather than model-driven. Our method is arguably simpler and easier to implement, and is equipped with rigorous analysis showing that it is an unbiased estimator of the underlying path. Practically, our method, like that of Cao and Krumm, also allows dierent directions of travel on a road to be identied, but since we do not use repelling forces, our method works for both left-hand and right-hand drive systems, and does not articially push lanes of opposite directions away from the centerline.

There have also been many developments in smoothing, noise-removal and outlier detection techniques in the GIS community: Nearest neighbor approaches for outlier detection are well studied [START_REF] Shekhar | Detecting graph-based spatial outliers[END_REF][START_REF] Shekhar | A unied approach to detecting spatial outliers[END_REF][START_REF] Lu | Multivariate spatial outlier detection[END_REF][START_REF] Chen | On detecting spatial outliers[END_REF][START_REF] Chen | Gls-sod: a generalized local statistical approach for spatial outlier detection[END_REF]. Chen et al. [START_REF] Chen | Leveraging spatio-temporal redundancy for rd data cleansing[END_REF] leverage spatio-temporal redundancy to cleanse RFID data. Guilbert

and Lin [START_REF] Guilbert | B-spline curve smoothing under position constraints for line generalisation[END_REF] developed a B-spline smoothing method with positional constraints for marine charts, where spatial conicts induced by normal B-spline smoothing

RR n°7754

can present navigational safety issues. Liu et al. [START_REF] Liu | Spatial outlier detection: random walk based approaches[END_REF], explore various random walk based approaches for spatial outlier detection. Arge, et. al. [START_REF] Arge | Cleaning massive sonar point clouds[END_REF], presented a method to clean sonar point clouds. They show theoretical analysis that their method is I/O ecient, and hence is practical for massive data sets.

Other related work in the GIS community includes that of Hönle et al. [START_REF] Hönle | Usability analysis of compression algorithms for position data streams[END_REF],

who presented an experimental study on various compression algorithms for trajectories, of Panangadan and Talukder [START_REF] Panangadan | A variant of particle ltering using historic datasets for tracking complex geospatial phenomena[END_REF], who use a data driven approach for tracking applications, of Johanssan and Jern [START_REF] Johansson | Geoanalytics visual inquiry and ltering tools in parallel coordinates plots[END_REF], who integrate statistical ltering techniques into a geographic visualization tool, and of Hall et al. [START_REF] Hall | Improving the quality of gps-based personal gazetteers[END_REF],

who investigate various heuristics to clean up personal GPS data.

In the broader research community, related work includes principal curves [START_REF] Hastie | Principal curves[END_REF], which are smooth curves passing through point clouds. They can be considered to be a nonlinear generalization of principal components. Cheng et al. [START_REF] Cheng | Cleaning uncertain data with quality guarantees[END_REF] use a data-driven approach to clean entries in a probabilistic database.

Chazal et al. [START_REF] Chazal | Geometric Inference for Probability Measures[END_REF] show connections between the distance to k-nearest neighbors and a natural distance function between probability measures. Their results can be used for robust topological and geometric reconstruction from data with heavy, but known, noise distributions. There has also been work on computing mean shapes [START_REF] Sánchez | A mean string algorithm to compute the average among a set of 2d shapes[END_REF][START_REF] Kendall | Shape Manifolds, Procrustean Metrics, and Complex Projective Spaces[END_REF] which study collections of shapes by considering their distributions in a shape-space.

Preliminaries

Recall that we have an underlying path, a database of traces and a query trajectory. We model the underlying path as a curve f : R → R D . We cannot observe this underlying path; instead, what we observe is the database of traces, which we model as noisy, discrete and nite vectors {f (i) } i with the following properties:

1. Each vector is of the form f (i) := (f

(i) αi , f (i) αi+1 , f (i) αi+2 , . . . , f (i) αi+m ) where f (i) x ∈ R D and f (i) x = f (x) + (i)
x where (i) x corresponds to an error term.

Each (i)

x is an independent and identically distributed error term with mean 0 and with E

(i) x 2 = σ 2 . Here, E (i)
x 2 is the trace of the covari- ance matrix of (i)

x , which characterizes the variances of (i)

x . In addition, we assume that there exists an upper bound M on the magnitude of the error, e.g.,

≤ M .

The query trajectory that the algorithm seeks to smooth is of the form g := (g αi , g αi+1 , g αi+2 , . . . , g αi+l ) where g x ∈ R D and g x = f (x) + η x where η (i) corresponds to an error term that has the same distribution (i)

x . The query trajectory may or may not come from the collection {f (i) } i . The smoothing objective is then to reduce the variance of the deviation errors of g with respect to the underlying path f by using the collection {f (i) } i , and output the smoothed trajectory as g.

Algorithmic Framework

The input to the algorithm consists of the traces {f (i) } i and a query trajectory g sampling an underlying path f . The algorithmic framework consists of three main steps:

1. Delay Embedding (Lifting): For the delay embedding, we choose an integer parameter n ≥ 0 that controls the length of the embedding. In a sense, this parameter controls the importance we place on the sequential information over the positional information. Now, given a vector where h x-m (i) is dened to be the ith coordinate of h x-m ∈ R D . For example, a vector ((1, 1), (2, 2, ), . . . , (k, k)) would be lifted to (

, . . . , (k -2, k -2, k -1, k -1, k, k)) for n = 1. (1, 1, 2, 2, 3, 3), (2, 2, 3, 3, 4, 4)
One caveat is the lifted vectors are now shorter, but in practice this is not much of an issue, since n can be set to value that is much smaller than the length of the vectors. We apply this lifting procedure to all points on f (i) and g to obtain f (i) and ĝ.

The advantages of using this lifting procedure are twofold. First, the sequential information around the point of interest encodes important attributes of the path such as speed and direction. This allows opposite directions of travel to be dierentiated, for example. Second, the delay embedding also allows us to take advantage of statistical concentration to reduce noise even further. In Section 4, we analyze this eect in detail.

2. Moving Towards Nearest Neighbors: We embed each point of the traces and query trajectory into the lifted space, forming a point cloud P in R (2n+1)D . Inspired by [START_REF] Chazal | Geometric Inference for Probability Measures[END_REF], in which points are moved along the gradient of the so-called distance-to-measure function, the smoothing technique works in the lifted space R (2n+1)D by moving the point ĝx to its nearest neighbor or to the barycenter of its k nearest neighbors in the point cloud P, i.e.,

gx = 1 k x∈knn(ĝx)
x where x, ĝx ∈ R (2n+1)D . . . . , gx+n (1), gx+n (2), . . . , gx+n (D))

RR n°7754

There are several alternatives for this procedure. One simple strategy is to use the center of gx , which is (g x (1), gx (2), . . . , gx (D)) as the coordinates for gx . We note that (g x (1), gx (2), . . . , gx (D)) may appear in dierent gx 's, which are not necessarily the same. Thus another strategy is that we can take the average of all (g x (1), gx (2), . . . , gx (D)) that appears in dierent gx 's as the coordinates for gx .

Our method is easy to summarize: we simply move the points of the lifted trajectory ĝx to the barycenter of several of its nearest neighbors in the lifted points of {f (i) } i . There are several variations of this framework: the particular one we will analyze is where we only look for one nearest neighbor from each trace f (i) . We reproject each lifted point to a point in the original space by averaging over all values of the smoothed coordinates of its image in the original embedding. For ease of reference, we outline the algorithmic framework in Algorithm 1.

Algorithm 1 Smoothing Framework

Input: Collection of traces {f (i) } i and a query trajectory g sampling an underlying path f . Output: Smoothed trajectory g. 1: Embed {f (i) } i and g in a high dimensional space as in Section 3 to obtain { f (i) } i and ĝ.

2: Move each point in ĝ towards several of its nearest neighbors in the points of {f (i) } i .

3: Recover trajectory g in original dimension by averaging relevant coordinates in ĝ.

Analysis

Smoothing Curves Using Delay Coordinates

In the discrete setting of this problem, we have information about each trace in {f (i) } i at discrete time steps. In practice, as in the example of a collection of GPS traces, the time stamps for dierent traces will not necessarily be synchronized, i.e. α i -α j is an integer for all i, j. For ease of analysis, however, we rst consider a simple case where each f (i) is synchronized with the query trajectory g. The analysis for the general asynchronized case will be provided in Theorem 4.5.

In the synchronized case, we can dramatically smooth the query trajectory (see Lemma 4.2, Theorem 4.3 and Proposition 4.4). This result is derived using the Hoeding inequality as described in Theorem 4.1, which was proved rst by Cherno [START_REF] Cherno | A measure of asymptotic eciency of tests of a hypothesis based on the sum of observations[END_REF] and Okamoto [START_REF] Okamoto | Some inequalities relating to the partial sum of binomial probabilities[END_REF] for the special case of binomial random variables:

Theorem 4.1 (Hoeding Inequality). Let X 1 , • • • , X n be independent bounded random variables such that X i falls in the interval [a i , b i ] with probability one, and

S n = X 1 + • • • + X n .
Then for any t > 0 we have

P S n -ES n ≥ t ≤ exp -2t 2 / n i=1 (b i -a i ) 2
When we embed a specic point g 0 on g, the lifted image ĝ0 is (g

-n , • • • , g n ).
We show that conditioned on f (i) being dened on indices -n to n, with high probability, the lifted image on f (i) whose delay coordinates have the same indices is the closest one. We have the following lemma: RR n°7754 Lemma 4.2. Assume that there exists a constant S ≥ 64σ 2 such that the underlying path

f k satises n k=-n f (k) -f (k + d) 2 ≥ Sd 2 n for every d ≥ 1. Let f (i) = f + (i) satisfy (i) ≤ M and E (i) 2 = σ 2 . Then, if n ≥ M 4 /σ 4 , with probability at least 1 -8 exp(-nσ 4 /M 4 ), the nearest lifted image on f (i) to (g n , g -n+1 , • • • , g n ) is (f (i) -n , f (i) -n+1 , • • • , f (i) n ).
Proof. For a piece of trace f (i) whose dening interval includes [-n, n], we consider the probability that (g

-n , • • • , g n ) is closest to (f (i) -n , f (i) -n+1 , • • • , f (i) n ). For xed d, we have P n k=-n g k -f (i) k 2 ≥ n k=-n g k -f (i) k+d 2 = P n k=-n η k - (i) k 2 ≥ n k=-n f k + η k -f k+d - (i) k+d 2 = P n k=-n 2(η k - (i) k+d ) T (f k+d -f k ) + n k=-n ( (i) k+d - (i) k ) T (2η k - (i) k - (i) k+d ) ≥ n k=-n f k -f k+d 2 ≤ P n k=-n (η k - (i) k+d ) T (f k+d -f k ) ≥ 1 4 n k=-n f k -f k+d 2 + P n k=-n ( (i) k+d - (i) k ) T (2η k - (i) k - (i) k+d ) ≥ 1 2 n k=-n f k -f k+d 2 Since (η k - (i)
k+d ) T (f k+d -f k ) are independent random variables that are bounded by 2M f k+d -f k , by concentration inequality, we know that

P n k=-n 2(η k - (i) k+d ) T (f k+d -f k ) ≥ 1 4 n k=-n f k -f k+d 2 ≤ exp - ( n k=-n f k -f k+d 2 ) 2 n k=-n 64M 2 f k -f k+d 2 RR n°7754
where the right hand side can be simplied as

exp - 1 64M 2 n k=-n f k -f k+d 2 Note that P n k=-n ( (i) k+d - (i) k ) T (2η k - (i) k - (i) k+d ) ≥ 1 2 n k=-n f k -f k+d 2 ≤ P n k=-n 2 η k 2 + 2 k 2 ≥ 1 2 n k=-n (f k -f k+d ) 2
Since Eη 2 k = σ 2 , then by concentration inequality again, we know the right hand side can be bounded by exp -

( n k=-n f k -f k+d 2 -8(2n + 1)σ 2 ) 2 16(2n + 1)M 4
Under the assumption that there exist a constant S such that

n k=-n f k - f k+d 2 ≥ Sd 2 n where S ≥ 36σ 2 , we have exp - 1 64M 2 n k=-n f k -f k+d 2 ≤ exp - Sn 64M 2 d 2
and by using the fact that S ≥ 64σ 2 we have exp -

( n k=-n f k -f k+d 2 -8(2n + 1)σ 2 ) 2 16(2n + 1)M 4 ≤ exp - 4(2n + 1)σ 4 M 4 (2d 2 -1) 2 Thus, P n k=-n (g k -f (i) k ) 2 ≥ n k=-n (g k -f (i) k+d ) 2 ≤ exp - Sn 64M 2 d 2 + exp - 4(2n + 1)σ 4 M 4 (2d 2 -1) 2 ≤ exp - nσ 2 M 2 d + exp - nσ 4 M 4 d
By summing up with respect to all d's, we know that the probability that (f

(i) -n , f (i) -n+1 , • • • , f (i) n ) is closest to (g -n , • • • , g n ) would be at least 1 - d =0 exp - nσ 2 M 2 d - d =0 exp - nσ 4 M 4 d

RR n°7754

Provided that n ≥ M 4 /σ 4 , the above quantity would be at least 1 -8 exp(-nσ 4 /M 4 ). Thus we conclude that the failure probability that the nearest lifted image on f (i) is not (f

(i) -n , f (i) -n+1 , • • • , f (i) n ) decays exponentially fast
with the number of delayed coordinates we use.

To interpret the conditions of the lemma, we note that n k=-n f (k)-f (k + d) 2 ≥ Sd 2 n requires that the path is always moving in some general direction. Note that if each component of f is strictly monotonic in some direction r with f •r always greater than a constant, then it is clear that we have

n k=-n f (k)- f (k + d) 2 ≥ Sd 2 n for all d ≥ 1 and n ≥ 1.
However, this is a much stronger condition than necessary, as f satisfying the condition can also have a kink or a zigzag pattern, provided that the size of the kink is not too large. Moreover, we require that S ≥ 64σ 2 , which means that the moving speed should not be too slow when compared with the level of noise. Lemma 4.2 bounds the probability that the delay coordinates of the nearest lifted image f (i) to (g -n , g -n+1 , • • • , g n ) are not centered around index 0 for any i. Thus, if a collection of N traces are given, we can use the union bound to bound the probability of identifying a wrong lifted image whose delay coordinates are not centered around 0, which leads to the following theorem: Theorem 4.3. Let {f (i) } i be a database of N noisy traces {f (i) } i , each following the assumptions of Lemma 4.2. Then with probability at least 1 -8N exp(-nσ 4 /M 4 ), the nearest lifted image on each

f (i) to (g n , g -n+1 , • • • , g n ) is (f (i) -n , f (i) -n+1 , • • • , f (i) n ), ∀i.
With the above theorem, it is easy to see that with high probability, we would use

1 N i f (i) 0
as the smoothed coordinate for g 0 . We now analyze the use of

1 N i f (i) 0
as an estimator for f 0 . Recall that if θ is an estimator of θ, then the bias and variance of θ is dened to be

Bias( θ) = E θ -θ, Var( θ) = E( θ -E θ) 2 .
Conditioned on the event A that the nearest lifted image on f (i) is identied for each i, which holds with probability at least 1 -8N exp(-nσ 4 /M 4 ), the conditional bias of the estimator

1 N i f (i) 0 is zero and the conditional variance of this estimator is σ 2 /N because E( 1
N i f (i) 0 |A) = f 0 , Var( 1 N i f (i) 0 |A) = σ 2 /N.
In summary, the algorithm has the following property in the synchronized case.

Proposition 4.4. With high probability, the smoothed coordinate of g 0 is an estimator for f 0 with zero bias. Furthermore, the variance of such an estimator can be reduced by a factor of N if we smooth the query trajectory g with respect to N relevant traces.

In the above theoretic analysis of the algorithm, we assume that each trace f (i) is synchronized with the query trajectory g. Proof. See Appendix A.

As we mentioned before, n k=-n f (k)-f (k +d) 2 is a stable measure of the moving speed of the path. In the case where S 1 and S 2 are far apart, one can still show that the lifted image on f (i) , whose delay coordinates are centered at the index k closest to 0, is one of the l-nearest neighbors of g -n , • • • , g n , where l depends on the dierence between S 1 and S 2 . In the asynchronized case, our estimator is not unbiased, but still succeeds in reducing the variance.

Connection with Moving Average

Our approach can be viewed of as a generalization of standard moving average techniques for smoothing trajectory data. Consider the case where a database of traces does not exist, and we move each point in ĝ towards several of its nearest neighbors as described in the second step of the smoothing framework. It is easy to see that we are essentially performing a moving average operation where we average the nearby coordinates of g to obtain a smoothed trajectory. More generally, we can assign weights on dierent dimensions in the lifted space, i.e., assigning weights w(r) on the rth dimension, we smooth g t as n r=-n w(r)g t+r .

We can further interpret the result of the moving average when a carefully choosen weighted kernel w is used. For example, if we choose an exponentially weighted kernel w(r) = √ λ 2 e - √ λ|r| , then it is possible to interpret the moving average result as an optimal solution which balances the geometry of the query trajectory g and the velocity constant moving trend in a continuous sense. More precisely, the moving average answers the following question:

Question 4.6. Given a noisy trajectory g : R → R. Can we generate a new trajectory y : R → R such that at each time t ∈ R, y(t) is close to g(t), while the speed y (t) is also close to a constant S?

Here, we focus on the one dimensional case, but it is straightforward to generalize the analysis to higher dimensions because dierent dimensions are smoothed independently in the moving average smoothing technique. We formulate this question as an optimization problem on a continuous nite interval [0, 1]:

min y 1 0 (y (t) -S) 2 dt + λ 1 0 (y(t) -g(t)) 2 dt (1)
If boundary conditions on y(0) and y(1) are given and we restrict the searching range of the trajectory y to the class of C 2 continuous functions on [0, 1], then using a calculus on variations approach we have the following lemma.

RR n°7754

Lemma 4.7. The optimal C 2 continuous trajectory y for the optimization problem (1) is characterized by the dierential equation

y (t) -λy(t) + λg(t) = 0 (2)
with given boundary conditions on y(0) and y(1).

Proof. See Appendix B.

We can solve the dierential equation (2) exactly using Laplace transforms: Lemma 4.8. All the solutions to the dierential equation (2) have the form

y(t) = c 1 e - √ λt + c 2 e √ λt + R √ λ 2 e - √ λ|r-t| g(r)dr (3)
In particular,

y(t) = R √ λ 2 e - √ λ|r-t| g(r)dr (4)
is the unique solution to (2) with the particular boundary values

y(0) = R √ λ 2 e - √ λ|r| g(r)dr y(1) = R √ λ 2 e - √
λ|r-1| g(r)dr.

Proof. See Appendix C.

When g is a noisy sample of the underlying path f , i.e., g(t) = f (t) + (t), the moving average technique generally yields a biased estimator [START_REF] Ruppert | Multivariate locally weighted least squares regression[END_REF] except for the case where the underlying trajectory f is linear. Specically, when f is linear, and is a mean zero continuous process, e.g., the generalized derivative of a Wiener process, then using y(t) =

= E R √ λ 2 e - √ λ|r-t| g(r)dr -f (t) 2 = E R √ λ 2 e - √ λ|r-t| w(r)dr 2 = R √ λ 2 e - √ λ|r-t| 2 dr = √ λ 4

RR n°7754

where the second to last equation follows from the properties of a Wiener process [START_REF]): an introduction with applications[END_REF]. Thus, when λ is large, we are essentially only using a few nearby coordinates around to smooth g; such a smoothing technique still has large variance which comes from the noises within g itself. When λ is small, however, such a smoothing technique works quite well.

As a nal remark, we note that the moving average can also be thought of as smoothing the query trajectory with respect to a collection of traces that consists of one single linear trace. To see this, consider the continuous formulation of the Algorithm 1 in the one dimensional case: for each point t ∈ R, we consider a small piece of g around t, e.g., the piece of g on [t -T /2, t + T /2]. The algorithm can then be interpreted as nding the piece of the standard trajectory f closest to this piece of g in the L 2 sense. Therefore we formulate a continuous version of the algorithm as the following:

min s T /2 -T /2 (g(t + r) -f (s + r)) 2 dr (5)
or, more generally, we can consider a weighted L 2 measure:

min s T /2 -T /2 w(r)(g(t + r) -f (s + r)) 2 dr (6)
If the standard trajectory f is linear, then without loss of generality we can assume f (s) = s. The result of the above optimization problem is exactly the moving average because the optimal solution s for the problem (6) satises d ds

T /2 -T /2
w(r)(g(t + r) -(s + r)) 2 dr = 0, which implies that the output of the smoothed g(t) is T /2 -T /2 w(r)g(t + r)dr T /2 -T /2 w(r)dr [START_REF] Bruntrup | Incremental map generation with GPS traces[END_REF] which is a moving average operation on g, as long as the weights w(r) are choosen to be a symmetric kernel, i.e., T /2 -T /2 w(r)rdr = 0.

Variations and Heuristics

Our smoothing framework naturally leads to algorithms for various smoothing objectives. One particular application we are interested in is to smooth an entire collection of GPS traces within a certain area as a noise-removal preprocessing step for algorithms that extract further information from the traces. To do this, we run our algorithm once for each trace in the collection, each time picking that trace as the query trajectory, and hence smooth all of traces in the collection. We use the small example of noisy traces shown Figure 2 to illustrate the eect of this application. Note that in this particular example, we are able to detect inaccuracies in the map, such as the shape of the upper-left cloverleaf from the collection of GPS traces. Unlike in the analysis and in the synthetic example, the GPS traces are usually noisy samples, not of a single underlying path, but of many dierent underlying paths derived by the road map and local trac laws. If we blindly apply the original algorithm, the result would be curves that appear smooth, but because traces far away from the query trajectory also aect it, the shape of the curves are sometimes distorted to a degree where the curve no longer follows the road network, see Figure 3. Note that in particular, near the bottom of the gure, one of the traces is moved to a region that is roughly in the middle of two roads. Theorem 5.1. Let ≤ M and η (i) ≤ M . Assume that the query trajectory g is sampled from f with g k = f k + k and h (i) = h + η (i) is another trace in the database where the lifted image f and ĥ have distance at least D, i.e., n k=-n f k -h k 2 ≥ D. Then, with probability at least 1 -exp(-D/128M 2 ), the distance between the lifted image ĥ(i) and ĝ is greater than D/2. Proof. See Appendix D. Thus, if the lifted images of two points from two dierent trajectories are suciently far apart, the thresholding ensures that their samples will not aect each other's smoothing results with high probability. Note that it is often possible to achieve the necessary separation between two road sections in the lifted space by simply increasing the parameter n.

RR n°7754

RR n°7754

Note that with the smoothed traces, it becomes much easier to distinguish dierent lanes and directions of travel than from the raw data, while the general shape of the curves is much better preserved than in Figure 3, in particular in the elongated cloverleaf pattern near the bottom of the page. This eect of the thresholding process is consistent throughout the entire data set.

Reparameterizing Data

We have noticed that the smoothing algorithm sometimes forms tiny loops or cusps these are a result of a vehicle slowing down dramatically and hence resulting in a point in the lifted space that is much further than it appears in

RR n°7754

Figure 5: Eect of thresholded algorithm after reparameterizing data.

the original space. Such features can be used to detect such anomalies in vehicle speed and direction that would have been less apparent in the original space.

If we are not interested in dierentiating traces with dierent speed as in the case of road network generation, we can reparameterize the input curves for constant speed, see Figure 5. Here, we see that the smoothed traces appear to have even less variance than the result in Figure 4. This can be attributed to the fact that without reparameterization, traces that appear to be close may have very dierent speeds, and hence their delay embedded points fall outside the threshold for inclusion into the smoothing algorithm.

Iterative Heuristic

In the event that we are smoothing an entire collection of traces, we can apply an iterative heuristic in which after we perform the smoothing step as analyzed in our framework, we iterate over the smoothed trajectories to further reduce error. That is, after we smooth each trajectory in the collection of traces, we use the smoothed trajectories as the new collection of traces to smooth. The eect of this smoothing approach is discussed in Section 6.4.

Experiments and Evaluation

The smoothing step of the algorithms relies heavily on algorithms for the nearest neighbor (NN) problem: Given a set of P points, one would like to preprocess P into a data structure such that given query point q, one could eciently nd the closest point in P to q. This is a classical problem in computational geometry [START_REF] Shamos | Closest-point problems[END_REF][START_REF] Clarkson | A randomized algorithm for closest-point queries[END_REF], but in high dimensions, theoretical results that improve on brute force search remain elusive. However if query answers are allowed to be approximate, there exist ecient solutions when the underlying space is Euclidean [START_REF] Indyk | Approximate nearest neighbors: Towards removing the curse of dimensionality[END_REF][START_REF] Har-Peled | A replacement for voronoi diagrams of near linear size[END_REF][START_REF] Sabharwal | Nearest neighbors search using point location in balls with applications to approximate voronoi decompositions[END_REF][START_REF] Andoni | Near-optimal hashing algorithms for approximate nearest neighbor in high dimensions[END_REF]. For more general metric spaces, there are algorithms for nearest neighbor search with a logarithmic query time, albeit with a constant

RR n°7754

that depends exponentially on the intrinsic dimension of the data [START_REF] Beygelzimer | Cover trees for nearest neighbor[END_REF]. In practice, ecient C++ libraries exist, such as the ANN Library by David Mount and Sunil Arya [START_REF] Mount | Ann: A library for approximate nearest neighbor searching[END_REF], which we use in our implementations. Our experiments were done in the context of smoothing an entire collection of trajectories, as this allows us to see the eects of the iterative heuristic. When the embedding dimension is low, our algorithm is fast in practice: On a 2.33GHz Macbook Pro with 3GB of RAM, we can smooth 7145 traces from the Taxicab dataset with a total of 54510 points using a lifted space of dimension 9 in several seconds. The result of this is shown in Figures 10 and13. The smoothing of the synthetic data set also took several seconds. For the Moscow dataset, we used a much higher embedding dimension and number of nearest neighbors, so the smoothing procedure took several hours. When the embedding dimension is high, a brute force nearest neighbor searching technique in lieu of the ANN library will most likely perform better. We note however, that the algorithm is fairly robust to parameter selection, and we can achieve results that are similar to the ones shown in Figures 9 and 12 in much less time.

Data Sets

To test the ecacy of the algorithm, we used two data sets: a synthetic collection of noisy and partial samples of a circular arc, and two real world data sets: One of 1331 GPS traces tagged Moscow from OpenStreetMap [2], and another of 7145 taxicab traces in an unindentied major city.

Synthetic Data Set To generate the synthetic data set, we discretized the circular arc from θ = 0 to 3π on the unit circle into a sequence of 100 points.

Then, to generate a individual trace, we selected a interval of this sequence uniformly at random, and then added an independent Gaussian noise with mean 0 and standard deviation of 0.1 to each point in this sequence. We repeated this process 100 times to generate 100 noisy and partial samples of the circular arc.

This data set is depicted in Figure 6.

Moscow Data Set To evaluate the performance of the algorithm on real world data sets, we used a set of 1331 GPS traces downloaded from Open-StreetMap, as shown in Figure 8. This set of GPS traces was obtained by cutting a larger collection of GPS traces in Moscow to the window between 55.74 • and 55.76 • latitude and 37.58 • and 37.62 • longitude. Our visualization was implemented using the QGraphicsScene framework of Qt, which would not display our traces until we scaled them up. We chose to scale longitude by a factor of 3486 and latitude by a factor of 24855, which also resulted in a reasonable projection for the Moscow area. The following discussion in this section will be using these scaled units when referring to the Moscow data set.

We then reparameterize the traces by sampling the original trajectories so that consecutive sample points are one units away from one another. A closeup of the data set can be seen in Figure 8, where it is evident that the traces are very noisy and have large variations in sampling error and density.

Taxicab Data Set A third data set we used to evaluate the algorithm is a set of 7145 traces obtained by cutting a larger of collection of Taxicab GPS traces

RR n°7754

Figure 6: Synthetic data: The original trajectories are in red while the trajectories smoothed using the algorithm are in dark red.

RR n°7754

RR n°7754

from an unidentied major city to a window of height 0.01612 • latitude and width 0.01662 • longitude. The taxicab traces were of varying sampling density and we only retained the sub-traces where consecutive sample points are less than 5 seconds apart and of distance less than 1000 meters apart. The latitude and longitude were both scaled by a factor of 1000, and consecutive samples were then taken every 0.1 units to reparameterize the traces.

Parameter Selection and Evaluation

We evaluated the results using both visual inspection and average distance to the underlying route. We selected parameters to produce reasonable results for both methods of evaluation. For visual inspection, we compared the smoothed traces to the underlying road network and veried that roads are not signicantly distorted, pulled apart from each other, or contracted. To compute the average distance to the underlying route, we used the approximate Fréchet map matching algorithm in [START_REF] Chen | Approximate map matching with respect to the fréchet distance[END_REF]. This algorithm nds a (1 + )-approximate closest path in the map with respect to the Fréchet distance. For two given curves

π 1 : [0, p 1 ] → R d and π 2 : [0, p 2 ] → R d , the Fréchet distance is dened as d F (π 1 , π 2 ) := inf α:[0,1]→[0,p 1 ] β:[0,1]→[0,p 2 ] max t∈[0,1] π 1 (α(t)) -π 2 (β(t)) 2
Intuitively, one can think of the Fréchet distance as the shortest leash possible for a man and a dog to walk on two curves, allowing their speeds to vary, but never walking backwards. For the ground truth map, we used maps downloaded from

CloudMade [1] and we set = 0.1. We note that the map matching is performed on the original trajectories and show that are smoothed trajectories are even closer to the matched routes than the original ones, even though the algorithm knows nothing about the underlying map. We measure distance to matched routes by using the average distance of the projection to the matched route.

The parameters we selected for the synthetic data set were n = 21 and 728 nearest neighbors while we selected n = 64 and 4096 nearest neighbors for the Moscow data set and n = 9 and 495 nearest neighbors for the Taxicab data set.

We did not need to select a threshold for the synthetic data set, while we used a threshold of 64 units for the Moscow data set and 2 units for the Taxicab data set.

Results

Visually, our results can be see in Figures 6, 9 and 10. We note that in the synthetic data set, our smoothing method produces traces along a smooth circle without contracting or otherwise shifting the centerline of the unit circle. We also note that in the Moscow and Taxicab data sets, error is reduced without distorting the curvature of individual roads. For the Moscow data set, we can now dierentiate dierent directions of travel, but they are not articially shifted apart. The Taxicab data set had much shorter traces, along with higher noise, so this eect was not as apparent with the parameters we selected. The results of the numerical performance metrics can be seen in Table 1. We note that for the Moscow and Taxicab data sets, on average, the traces were already very close to the map, and given that the map we used may not be entirely accurate, our results present a reasonable reduction in error.

RR n°7754

RR n°7754

Figure 10: Smoothed Taxicab Data Set.

Iterative Heuristic

We also tested the iterative heuristic by running the smoothing algorithm through four iterations, using the same parameters as before. The results are shown in Table 1 and Figures 11,12, and 13. Interestingly, although the Taxicab data set appears much smoother after the iterative heuristic, the average distance to the matched route actually increased slightly. Possible reasons for this phenomenon include drift in the iterative smoothing process, inaccurate maps, or biased sampling from the original data set. Nevertheless, the iteratively smoothed collection of traces appears much cleaner, which can be useful for many applications.

Conclusion

We have presented a simple and practical algorithmic framework for smoothing a query trajectory with respect to a collection of traces. Under natural assumptions, our method can be proven to reduce variance among a noisy collection of traces. Our framework can also be viewed of as a generalization of the classical moving average technique, and we show additional connections to that method as well. Experimentally, we show that the algorithm works well on both synthetic data, and on two dierent collections of GPS traces. We also explore variations in the algorithm, and demonstrate how they result in dierent smoothing behavior. There are several directions for future work these

RR n°7754

Figure 11: Synthetic data after 4 smoothing iterations.

RR n°7754

RR n°7754

include inferring a reasonable number of delay coordinates and threshold value from the data itself without the use of a ground truth map, as well as examining special structure in the lifted images that might facilitate faster nearest neighbor queries.

Generally, we consider the probability that the lifted image g

-n , • • • , g n is closer to f (i) k-n , f (i) k-n+1 , • • • , f (i) k+n than f (i) k+d-n , • • • , f (i) k+d+n . Then we have P n j=-n g j -f (i) k+j 2 ≥ n j=-n g j -f (i) k+j+d 2 = P n j=-n f j + η j -f k+j - (i) k+j 2 ≥ n j=-n f j + η j -f k+j+d - (i) k+j+d 2 = P n j=-n 2(η j - (i) k+j ) T (f j -f k+j ) + n j=-n 2(η j - (i) k+j+d ) T (f k+j+d -f j ) + n j=-n ( (i) k+j+d - (i) k+j ) T (2η j - (i) k+j - (i) k+j+d ) ≥ n j=-n f j -f k+j+d 2 - n j=-n f j -f k+j 2 Denote n j=-n f j -f k+j+d 2 - n j=-n f j -f k+j 2
by t, then we have the above probability is bounded by

P n j=-n 2(η j - (i) k+j ) T (f j -f k+j ) ≥ t 3 + P n j=-n 2(η j - (i) k+j+d ) T (f k+j+d -f j ) ≥ t 3 + P n j=-n ( (i) k+j+d - (i) k+j ) T (2η j - (i) k+j - (i) k+j+d ) ≥ t 3 Since 2(η j - (i)
k+j )(f j -f k+j ) are mean zero independent random variables and they are bounded by 2M f j -f k+j , by concentration inequality, we know that

P n j=-n 2(η j - (i) k+j ) T (f j -f k+j ) ≥ t 3 ≤ exp - t 2 96 n j=-n M 2 (f j -f k+j ) 2 (8)

RR n°7754

For the same reason, we have

P n j=-n 2(η j - (i) k+j+d ) T (f k+j+d -f j ) ≥ t 3 ≤ exp - t 2 96 n j=-n M 2 (f j -f k+j+d ) 2 (9) Note that P n j=-n ( (i) k+j+d - (i) k+j ) T (2η j - (i) k+j - (i) k+j+d ) ≥ t 3 ≤ P n j=-n 2 η j 2 + 2 k+j 2 ≥ t 3 (10)
Now we consider the two cases:

1. The nearest time point k to 0 satises |k| < 1/4. In such a case, let |d| ≥ 1, then under the assumption that 2S 1 ≥ S 2 we can easily verify that

t = n j=-n f j -f k+j+d 2 - n j=-n f j -f k+j 2 ≥ S 1 (k + d) 2 -S 2 k 2 ≥ 1 2 S 1 d 2
Thus, using concentration inequality we know (8) and ( 9) are both bounded by exp(-C 1 S 1 d 2 n/M 2 ) where C 1 is a constant. [START_REF] Chen | Approximate map matching with respect to the fréchet distance[END_REF] can also be bounded by exp(-C 2 σ 4 dn/M 4 ) for the same reason as in the proof of Lemma 4.2. By summing up the tails, we know that with probability at least 1 -8 exp(-Cσ 4 n/M 4 ), the lifted image whose delay coordinates centered at x on f (i) is the nearest neighbor of g -n , • • • , g n . Proof. The general solution to the dierential equation ( 2) is given by y(t) = y p (t)+y h (t) where y p (t) is a particular solution of the nonhomogeneous equation RR n°7754

(2) and y h (t) is the general solution of the associated homogeneous equation y (t) -λy(t) = 0. Proof. We consider the probability that (g -n , • • • , g n ) and (h (i)

-n , • • • , h (i)
n ) is less than D:

P n k=-n g k -h (i) k 2 ≤ D/2 = P n k=-n η (i) k -k 2 + f k -h k 2 + 2( k -η (i) k )(f k -h k ) ≤ D/2 ≤ P n k=-n 2( k -η (i) k )(f k -h k ) ≤ D/2 -D ≤ exp -2(D/2 -D) 2 / n k=-n 64M 2 f k -h k 2 = exp -(D/2 -D) 2 /32M 2 D = exp -D/128M 2 RR n°7754

3 .

 3 Recovering the Trajectory: With each point ĝx moved to a new position (denoted by gx ) in the last step, we now recover a new trajectory g from the points ĝx 's. Let gx = (g x-n (1), gx-n (2), . . . , gx-n (D); gx-n+1 (1), gx-n+1 (2), . . . , gx-n+1 (D);

  g(r)dr as an estimator for f (t) has zero bias because E f (r)dr = f (t) is only true for the case where f is linear. Moreover, the variance for such an estimator is Var

Figure 2 :

 2 Figure 2: Noisy Traces.

Figure 3 :

 3 Figure 3: Eect of unmodied algorithm.

Figure 4 :

 4 Figure 4: Eect of thresholded algorithm.

Figure 7 :

 7 Figure 7: Moscow Data Set.

Figure 8 :

 8 Figure 8: Taxicab Data Set.

Figure 9 :

 9 Figure 9: Smoothed Moscow Data Set.

Figure 12 :

 12 Figure 12: Moscow data set after 4 smoothing iterations.

Figure 13 :

 13 Figure 13: Taxicab data set after 4 smoothing iterations.

2.f j -f k+j 2 ≥ S 1 S 1 d 2 1 0 1 0η 1 0 1 0

 2121111 The nearest time point k to 0 lies outside of [-1/4, 1/4]. In such a case, without loss of generality, we can assume k < 0.When d < -1 or d > 2, (k + d) 2 -S 2 k 2 ≥ 1 2Thus, for any d < -1 or d > 2, (8) and (9) are bounded by exp(-C 1 S 1 d 2 n/M 2 );[START_REF] Chen | Approximate map matching with respect to the fréchet distance[END_REF] is bounded by exp(-C 2 σ 4 dn/M 4 ) for certain constants C 1 and C 2 , by summing up of the tails, we know that with at least 1 -8 exp(-Cσ 4 n/M 4 ), the lifted image whose delay coordinates centered at x on f (i) will be one among the two nearest neighbor ofg -n , • • • , g n . t) -S) 2 dt + λ (y(t) -g(t)) 2 dtas A(y). Suppose η : [0, 1] → R is a C 2 function and vanishes at the end points 0 and 1. Then we must have:A(y) ≤ A(y + η)for any value close to 0. Therefore the derivative of A(y + η) with respect tomust vanish at t) + η (t) -S) 2 | =0 dt t) + η(t) -g(t)) 2 | =0 dt = (t)(y (t) -S)dt + λ η(t)(y(t) -g(t))dt = η(t)(y (t) -S))(λy(t) -λg(t) -y (t))dtThus, η(t)(λy(t) -λg(t) -y (t))dt = 0 for any twice dierentiable function η that vanishes at the endpoints. We can now apply the fundamental lemma of calculus of variations: If b a η(t)H(t)dt = 0 for any suciently dierentiable function η(t) within the integration range that vanishes at the endpoints of the interval, then it follows that H(t) is identically zero on its domain. Hence λy(t) -λg(t) -y (t) = 0. C Proof of Lemma 4.8

  The solution to the homogeneous equation is given by c1 e - √ λt + c 2 e √ λt .Now we verify that y p (t) = R √ λ 2 e - √ λ|r-t| g(r)dr is a particular solution to (2).To show this, we note that y p (t) = r)dr -λg(t).

  The last identity is equivalent to y p (t) = λy(t) -λg(t).

  h = (h α , h α+1 , . . . h α+m ), we embed each h x ∈ R D into the lifted space R D(2n+1) by mapping it to its lifted image ĥx :

	ĥx = (h x-n (1), h x-n (2), . . . , h x-n (D);
	h x-n+1 (1), h x-n+1 (2), . . . , h x-n+1 (D);
	. . . , h x+n (1), h x+n (2), . . . , h x+n (D))

  Theorem 4.5. Assume that there exist constants S 1 , S 2 and B such that S 1 ≥ Bσ 2 for which the pathf k satises -f (k + d) 2 ≤ S 2 d 2 nfor every d ≥ 1, where S 1 and S 2 further satises 2S 1 ≥ S 2 . If n ≥ M 4 /σ 4 , then there exists a constant C which depends only on B, such that with probability at least 1 -8 exp(-Cnσ 4 /M 4 ), the lifted image on f (i) centered at the index k closest to 0 is one of the two nearest neighbors of (g -n , • • • , g n ).

	n
	S 1 d 2 n ≤	f (k)
	k=-n

However, in real applications, this is rarely true. By carrying out a similar analysis as in Lemma 4.2, we can show the following: RR n°7754

Table 1 :

 1 Fréchet distance to routes obtained using map matching on the original trajectories. The distances for the Moscow and Taxicab data sets are using longitude/latitude coordinates. Approximate upper bounds in meters for the Moscow data set are shown.

		Synthetic	Moscow	Taxicab
	Mean Dist	0.0795	0.000143	0.000196
	Orig		< 15.9m
	Mean Dist	0.00553	0.000114	0.000156
	Smoothed		< 12.7m
	Mean Dist	0.00461	0.000111	0.000174
	Iter Smoothed		< 12.3m

RR n°7754

http://www.inria.fr ISSN 0249-6399

Acknowledgments

This work was supported in part by ARO grant W911NF-10-1-0037, the AH-PCRC Center at Stanford (ARO grant W911NF-07-2-0027), NSF CCF grant 1011228, a Google Faculty Research grant, GIGA grant ANR-09-BLAN-0331-01, and the European project CG-Learning No. 255827. We would also like to thank the associated research team COMET.


Di Fiore

Tom Van

Laerhoven Fabian
email: tom.vanlaerhoven@uhasselt.be

William Van
email: william.vanhaevre@uhasselt.be

Haevre Frank

Van Reeth

Computer Animation and Virtual Worlds

Keywords: animation, 2D animation, painterly animation, paint system, painton-glass animation, 2D animation, painterly animation, paint system, paint-on-glass

In this paper we deal with paint-on-glass animation, which is a technique for mak-1 Page 1 of 25 http://mc.manuscriptcentral.com/cavw -For Peer Review Computer Animation and Virtual Worlds ing animated films by pushing slow-drying paints from frame to frame directly under the camera. As artwork is continuously destroyed upon creating new frames, the animator is not able to rehearse or refine the animation afterwards. Furthermore, due to impracticable issues like how to stack up layers containing wet paint or how to overlay masks on the wet medium, one has to take both creative as technical decisions for each shot.

Our approach consists of an interactive paint setup that physically simulates paint media. Together with a set of digital tools the artist is relieved from the difficult task of sustaining a constant frame-to-frame coherence while animating and is given the possibility to modify or undo earlier paint modifications. Regarding the interaction part, the artist stays in full control by employing a tangible interface. This allows artists to push pigment with their fingers, use special-purpose digital brushes, or even to employ real life tools like paper tissue.

Feedback from artists confirms the real-life behaviour of stretching and smudging the paint as well as interacting with the setup, resulting in a natural and reality-based methodology.

Introduction

In spite of its name, paint-on-glass animation is best defined as a technique for making animated films by pushing around some kind of (wet) medium directly under the camera and recording it frame-by-frame (Figure 1). It is categorised as an "under-the-camera-technique" as it is nearly always undertaken by an individual artist or animator rather than by a production studio using factory-like processes. As the medium is pushed around directly under the camera and recorded frame-by-frame, each image seems to merge from the previous one and melt into the next, resulting in movement that can be very fluid and organic -a continual process of metamorphosis.

A well-known practitioner of the paint-on-glass animation technique is Russian animator Aleksander Petrov. His most famous work is the award winning short film "The Old Man and the Sea" which took two and a half years to complete [START_REF] Petrov | World Wide Web[END_REF] illustrating the complexity of the technique.

Motivation. As the medium is pushed around from frame to frame, artwork is continuously destroyed as new artwork is created. Without the ability to rehearse and refine the animation and without the disposal of layering and masking tools, the animator must plough on regardless incorporating any errors into the sequence.

On the one hand this stimulates spontaneity and delivers work that is very fresh and distinctive as it celebrates the method of its making. For example, characters may move from place to place not by walking, but by being smudged away to re-form out of the background 3 at the required location. On the other hand, this particular form of step-ahead animation is a painstaking process which demands much self-discipline, endurance and concentration as the resulting frames are never a series of clean images of the kind produced by other animation techniques, but also contain a record of their making due to the pushing process.

Contribution.

In this paper we present our approach to facilitate paint-on-glass animation. In close collaboration with artists and animators we developed a set of digital tools that do enable them to achieve the fine quality of a finished painting for each frame, while preserving all artisan control. Primarily, our system features:

• a real-time free-form paint tool to physically simulate the complex behaviour of different paint media (gouache, watercolour, impasto, pastel);

• painting/animating by pushing the drawing media around resulting in a stretch and smear effect;

• a tangible interface consisting of a multitouch screen allowing the use of one's fingers as well as physical media including brushes and tissues;

• several tools that are absent in the traditional method due to the wet nature of the painting media such as working with layers and masks;

• the availability of animation tools like onion skinning and flipbook animation; Paper Organisation. First, we survey work we consider related to ours. Next, we describe our approach starting from the employed painting setup to an overview of the digital painton-glass animation process. Finally, clarifying examples are provided in which we also discuss our results.

Related Work

We start this section by giving a brief overview of traditional paint-on-glass animation.

Then, we elaborate on digital paint solutions that we consider related to ours.

Traditional Paint-on-Glass Animation

What we give here is an account as it might be for a practiced artist in a home studio.

We separately discuss the painting process (i.e. painting a single frame) and the animation process (i.e. making a series of drawings).

Painting. The painting process itself resembles painting on a canvas as the initial paint is added with brushes or other tools including sponges, cotton buds, small sticks, strong tissue or even fingers. As the paint medium will be reused from frame to frame, the artist is constrained to work with slow drying paint media like oil paint in order to keep working with the medium across several hours or even days. This causes also many impracticable issues like how to stack up layers of wet paint or how to overlay masks on the wet medium.

Animating. For classic 2D animation the animator draws or sets up objects one frame at a time in sequential order until the sequence is complete, also referred to as step-ahead animation. In this way there is one physical drawing or image per frame. In the case of paint-on-glass animation the artist is working directly under-the-camera: the artist has to record the current frame before proceeding to the next one as current artwork is continuously destroyed for creating each new frame. When a sequence has been filmed, there is no raw material left besides the final film. So, while animating one cannot rely on techniques like onion skinning or flipbook animation as previous artwork is not available anymore. Hence, there is also no artwork to turn to when modifying the film afterwards.

Digital Paint Systems

In this section we will discuss digital paint solutions that we consider related to ours.

Painting. Looking at the majority of published paint systems [START_REF] Saito | 3D physics-based brush model for painting[END_REF][START_REF] Baxter | DAB: Interactive haptic painting with 3D virtual brushes[END_REF][START_REF] Chu | Real-time painting with an expressive virtual Chinese brush[END_REF][START_REF] Van Laerhoven | Brush up your painting skills: Realistic brush design for interactive painting applications[END_REF], virtual brushes are by default operated by some kind of digital stylus in combination with a digital tablet.

Although popular and basically intuitive, this way of working is nothing more than an ap- Other paint systems explicitly target novel interfaces to bridge the gap between physical and digital painting. Vandoren et al. developed IntuPaint, which is a tangible interface for a digital paint easel, using an interactive surface and electronic brushes with a tuft of bristles [START_REF] Vandoren | IntuPaint: Bridging the gap between physical and digital painting[END_REF]. This approach provides natural interaction and enables detailed tracking of specific brush strokes. Additional tangible and finger-based input techniques allow for specific paint operations or effects. Okaichi et al. propose an interactive oil paint system using a painting knife [START_REF] Okaichi | A virtual painting knife[END_REF] that is controlled using a haptic feedback device. Microsoft Surface [START_REF]Microsoft surface. World Wide Web[END_REF] is based on the IR image capture of an IR illuminated diffuser screen underneath a transparent surface allowing natural finger-based interaction. Presence of real objects (e.g., brushes) can be detected, but only a blurred image of the brush is generated. Recently, Vandoren et al.

presented a novel interactive canvas FluidPaint allowing the use of real wet brushes [START_REF] Vandoren | Review Computer Animation and Virtual Worlds painting system using real wet brushes[END_REF].

Animating. Regarding animation, only small initiatives have been taken so far. Van Laerhoven et al. introduce the concept of brushes and pigments enhanced with behavioural intelligence that allow the user to enrich and animate interactively created images by semiautonomously embedding procedural animation into them [START_REF] Van Laerhoven | Hand-painted animation with intelligent brushes[END_REF]. Some commercial digital art applications also claim to target at creating animations ( [START_REF] Painter | World Wide Web[END_REF][START_REF]Project Dogwaffle[END_REF][START_REF]Adobe PhotoShop[END_REF]), however, the animation part mainly is limited to the availability of a keyframe editor and the use of animated brushes to draw frames along a path.

Approach

In this section we give an overview of the main parts of our approach. We first describe the employed digital paint setup. Next, similar to the traditional process of previous section we will highlight both the painting and animation process.

Digital Paint Setup

For the simulation part, we make use of an interactive paint system that adopts physicallybased algorithms to simulate the complex interaction of a layered virtual canvas with different paint media, including watercolour, gouache and impasto, but also dry media like pencil or pastel [START_REF] Van Laerhoven | Real-time simulation of watery paint[END_REF]. This paint system encapsulates the drying process, paint diffusion as well as different pigment attributes (Figure 2(a)). Furthermore, it relies on Navier-Stokes for defining the velocities of the fluid body and on the specialised Kubelka-Munk colour mixing scheme for rendering paint layers. These features are vital elements to obtain the much sought after 'natural look' in digital artwork. Every step of the simulation makes full use of programmable graphics hardware to achieve interactive simulation rates.

Another important component of the system relates to its virtual brush design that relies on constrained energy optimisation and on an anisotropic friction model to deform the brush tuft and adequately translate the user's input (e.g., originating from a tablet interface) into realistic and predictable strokes. This design also allows for bidirectional pigment transfer between the canvas and the selected drawing tool. Regarding the interaction part, we chose to use a tangible interface to enable users to intuitively manipulate the drawing media (Figure 2(b)). The setup consists of an interactive multitouch surface based on the principle of frustrated total internal reflection (FTIR) [START_REF] Vandoren | IntuPaint: Bridging the gap between physical and digital painting[END_REF].

The contact surface of each object on the paint table is tracked by a 0.3 megapixel video camera, equipped with an infrared band-pass filter, after which the footprint of the objects touching the surface is transferred to the paint simulation for further processing and rendering. In contrast to the original paint table, our setup is able to track multiple objects on the surface simultaneously. This way, artists can push pigment with multiple fingers at once, use special-purpose digital brushes with tufts, or even employ paper tissue.

A major advantage of this setup is that it provides coinciding action and perception spaces, resulting in a natural, reality-based methodology. In addition, it encourages multiple users to work collaboratively on the same artwork.

Painting

Starting from a blank canvas, the artist creates or alters an image by adding some kind of paint medium (e.g., gouache, impasto, watercolour, pastel) using fingers or several available tools (brushes, pencils, spatulas, sponges, fingers). A background image can be depicted as a reference and parts of the drawing canvas can interactively be masked to prevent paint from reaching or settling down on the masked parts. Drawings can also be separated into layers allowing for robust manipulations of the frame content such as depth ordering and hid- Watercolour images made with the system contain the typical effects that can be recognised in images produced with real thin paint, like the dark-edge effect, watercolour glazing, wet-on-wet painting and the use of different pigment types.

Animating

While animating, each subsequent frame in the animation is created by pushing the drawing medium of the artwork in the current frame towards its desired position for the next frame.

During interaction, movement on the multitouch surface is converted to a vector field based on the current velocity of the brush/finger (see Figure 3(a)), which is then applied to the paint (both pigment and water) by means of an advection step -(⃗ ⋅ ∇) , where ⃗ denotes the 2D velocity vector field and the amount of paint at a certain position on the canvas. As in our paint system mass is contained in the grid cells of the digital canvas, the advection scheme boils down to individual cells exchanging content with their neighbours (Figure 3(b)) [START_REF] Van Laerhoven | Real-time simulation of watery paint[END_REF].

For each cell, we measure the volume of paint that is exchanged with all neighbouring cells as follows.

When calculating the amount of mass that flows from the centre cell to its right neigh-

⃗ average = ⃗ center + ⃗ right 2 (1)
Within a given time step Δ , a distance Δ = (⃗ average ) Δ is travelled. In case we are calculating the movement of water, the area covered by the volume of displaced water is given by the coloured area in Figure 3 water quantity is given by the following equation:

Δ = Δ (cellwidth × cellheight) ( 2
)
This procedure is repeated to calculate the exchanged water quantities with the three remaining neighbours. We add up the results and divide it by four, because each neighbour contributes exactly one quarter to the total flux. The same approach can be followed to work out the movement of pigment.

Figure 4 depicts some stills of a user pushing paint (in this case impasto) with two fingers on an interactive multitouch surface.

In order to assure more control over the way the artist is obliged to make things move, he/she can rely on digital animation tools to make manipulative operations more predictable 11 including a keyframe editor, onion skinning and a flipbook tool (Figure 5). The onion skinning tool aids the artist to achieve good frame-to-frame coherence, resulting in smoother frame transitions by showing trace layers of previously painted frames whereas the digital flipbook tool provides immediate feedback on the status of the animation being created.

Results

Figure 7 depicts some snapshots of a paint-on-glass animation of a bird created with impasto paint illustrating the stretch and smear effects. Even for the start frame, the artist chose to gradually refine a basic image (painted using a brush) using his fingers only (see Figure 6). The entire animation consists of 31 frames of which each is a modified version of its predecessor by only pushing paint from frame to frame.

Discussion. The target audience of our setup is skilled artists as the paint-on-glass technique itself requires some experience in the field of painting. The platform itself is targeted at both traditional and digital artists as the interface is similar to the traditional environment.

Therefore traditional as well as digital painters were involved in an informal user test around 20 frames a day, on average). Although the test session was rather short we can say that all participants were very positive. They found the user interface familiar and praised the real-life behaviour of manipulating the paint as well as interacting with the setup, resulting in a natural and reality-based methodology.

Conclusions

In this paper we dealt with making animated films by pushing around wet media directly under the camera, known as paint-on-glass animation. Without the ability to rehearse and refine the animation as opposed to classic 2D animation, this particular form of animation is a painstaking process for which animators need to persevere regardless incorporating any errors or unwanted smudging into the sequence.

Our approach tackles the issues inherent to paint-on-glass animation by offering the animator an interactive paint system that adopts physically-based algorithms to simulate the complex interaction of a virtual canvas with different paint media, and a set of digital tools that do enable him/her to achieve the fine quality of a finished painting for each frame while preserving all artisan control. Regarding the interaction part, a tangible interface was chosen to enable users to intuitively manipulate the drawing media using their own hands.

13

Future Work. The emphasis in this paper is on paint-on-glass animation. However, the same approach could be applied to any "under-the-camera-technique" such as animating using fine sand, sugar or coffee grind. For the simulation part this requires a bidirectional transfer of particles; particles need to be deposited or scraped off again. Regarding the interaction process we have to consider how to sprinkle dry media.

Furthermore, we strongly believe in investigating new interaction techniques which are not present at all in traditional animation such as using simple hand gestures and finger taps to easily navigate between frames or assets. Inspired by ShapeTouch [START_REF] Cao | ShapeTouch: Leveraging contact shape on interactive surfaces[END_REF] another option is to exploit the shape of hands to work on the interactive surface just like when dealing with real physical objects (e.g., pinning, peeling or flicking).

http://mc.manuscriptcentral.com/cavw -For Peer Review

4 Page 4

 44 of 25 http://mc.manuscriptcentral.com/cavw -For Peer Review Computer Animation and Virtual Worlds • digital tools to control the drying time, correct mistakes, remove pigment and/or water as desired, and the ability to save, load or reuse (intermediate) results.

5 http:

 5 //mc.manuscriptcentral.com/cavw -For Peer Review Computer Animation and Virtual Worlds

6 Page 6

 66 of 25 http://mc.manuscriptcentral.com/cavw -For Peer Review Computer Animation and Virtual Worlds proximation of 'the real thing'. Issues such as parallax, the lack of tufts and absence of physical tools need to be solved first.

7 Page 7

 77 of 25 http://mc.manuscriptcentral.com/cavw -For Peer Review Computer Animation and Virtual Worlds

8 http:

 8 //mc.manuscriptcentral.com/cavw -For Peer Review Computer Animation and Virtual Worlds

9 http:

 9 //mc.manuscriptcentral.com/cavw -For Peer Review Computer Animation and Virtual Worlds ing/revealing objects. Furthermore, users can interactively control aspects like drying time (including the possibility to keep paint permanent wet or pause the drying), undo mistakes, remove pigment and water as desired and have the ability to save, load or reuse intermediate results.

10 Page 10

 1010 of 25 http://mc.manuscriptcentral.com/cavw -For Peer Review Computer Animation and Virtual Worlds bour with velocities ⃗ center and ⃗ right respectively, as shown in Figure 3(b) (left), any point along this border has velocity:

  (b) (left), and equals Δ (cellheight). The total volume of displaced water then is Δ = Δ (cellheight) i,j (Figure 3(b) (right)). The change in

12 Page

 12 evaluating the setup as well as creating the examples. Only basic instructions were given (i.e. explaining the user interface). After that, the test participants were allowed to work on their subjects of choice without any training.All participants worked between two and three hours to complete their animations (43 mc.manuscriptcentral.com/cavw -For Peer Review Computer Animation and Virtual Worlds frames on average) which is significantly faster than the traditional way of working (i.e.

Figure 1 :Figure 2 Figure 3

 123 Figure 1: An artist at work using an "under-the-camera-technique" with a combination of

Figure 4 :

 4 Figure 4: Some stills depicting a user pushing paint (impasto) with two fingers on an inter-

Figure 5 :Figure 6 :

 56 Figure 5: When drawing a new frame the animator can check temporal coherence by flipping

Figure 7 :

 7 Figure 7: Some snapshots of a paint-on-glass animation of a bird created with impasto

Page 11 of 25 http

://mc.manuscriptcentral.com/cavw -

For Peer Review Computer Animation and Virtual Worlds

Acknowledgements

We gratefully express our gratitude to the European Fund for Regional Development and the Flemish Government, and to Xemi Morales and Karel Robert for their artistic contribution.