Denis Fortun
email: denis.fortun@epfl.ch

Patrick Bouthemy
email: patrick.bouthemy@inria.fr

Charles Kervrann
email: charles.kervrann@inria.fr

A Variational Aggregation Framework for Patch-Based Optical Flow Estimation

Keywords: Optical flow, parametric motion, aggregation, variational optimization

We propose a variational aggregation method for optical flow estimation. It consists of a two-step framework, first estimating a collection of parametric motion models to generate motion candidates, and then reconstructing a global dense motion field. The aggregation step is designed as a motion reconstruction problem from spatially-varying sets of motion candidates given by parametric motion models. Our method is designed to capture large displacements in a variational framework without requiring any coarse-to-fine strategy. We handle occlusion with a motion inpainting approach in the candidates computation step. By performing parametric motion estimation, we combine the robustness to noise of local parametric methods with the accuracy yielded by global regularization. We demonstrate the performance of our aggregation approach by comparing it to standard variational methods and a discrete aggregation approach on the Middlebury and MPI Sintel datasets.

Introduction

General positioning

Optical flow estimation is based on the conservation assumption of image features such as image intensity, image gradient or texture descriptor. The so-called brightness constancy assumption is the most used one. It provides a single equation and is consequently insufficient to recover the two components of the motion vector. A usual way to overcome this under-determination is to impose a spatial coherency constraint for the flow field. Existing methods can be classified into two main categories:

-Local spatial coherency is exploited when considering a parametric motion model, e.g., local translation [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF], affine model or quadratic model [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF], in a given neighborhood or an appropriate local region. The neighborhoods must be sufficiently textured or contain interest points to supply reliable velocity vectors. -Global coherency [START_REF] Horn | Determining optical flow[END_REF] imposes a regularization constraint to the motion field on the whole spatial domain. The flow field is generally assumed to be piecewise smooth and the strategy is to minimize a global energy of the form

E(w) = Ω ρ(x, u, v, w) + λ φ(∇w(x)) dx, (1)
that explicitly combines a potential ρ(•), which penalizes deviations from the brightness constancy equation, with a regularization potential φ(•) which penalizes high values of the norm of the gradient ∇w of the velocity field w : Ω → R 2 , where Ω ⊂ R 2 denotes the image domain. The two consecutive images are denoted by u, v : Ω → R, x ∈ Ω denotes one pixel of the grid Ω, and λ is a balance parameter.

The best state-of-the-art results are obtained with the global approach. Nevertheless, several open issues remain unsolved. One of the main limitations comes from the undesirable effects due to the coarse-to-fine strategy used to handle large displacements [START_REF] Enkelmann | Investigations of multigrid algorithms for the estimation of optical flow fields in image sequences[END_REF]. The motion of small objects is discarded at coarse scales, and the error is propagated in the incremental updates at finer scales when the displacement is larger than the object size. As a result, motion details are not correctly recovered in the final estimated flow field [START_REF] Xu | Motion detail preserving optical flow estimation[END_REF]. Large displacements are also associated to large occlusions, which are another major source of errors. Occlusion handling is often treated as a post-processing task. It is then very sensitive to errors in the initial motion estimation. Finally, noise sensitivity is usually ignored in standard optical flow evaluation benchmarks. However, if pixelwise data potentials provide best results in absence of noise, they are not adapted when noise is present in input images. To limit the impact of these failure cases, the solution of the global approach is often to increase regularization, producing oversmoothed results, loosing motion details and blurring discontinuities.

Existing purely local methods [START_REF] Bigun | Multidimensional orientation estimation with applications to texture analysis and optical flow[END_REF][START_REF] Jodoin | Optical-flow based on an edge-avoidance procedure[END_REF][START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF][START_REF] Senst | Robust local optical flow for feature tracking[END_REF] are far from being able to compete with global methods in terms of accuracy in optical flow benchmarks. The main issue is to be able to select appropriate local regions. The most basic approaches considering square patches centered on each pixel [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF] are unable to retrieve motion discontinuites. They are also prone to the same large displacement and occlusion problems as the global methods. Nevertheless, joint global motion estimation and segmentation approaches [START_REF] Mémin | Dense estimation and object-based segmentation of the optical flow with robust techniques[END_REF][START_REF] Sun | Layered segmentation and optical flow estimation over time[END_REF][START_REF] Unger | Joint motion estimation and segmentation of complex scenes with label costs and occlusion modeling[END_REF] have demonstrated that piecewise parametric representation of flow fields can yield excellent results when local regions are appropriately chosen. However, the required alternate optimization scheme is computationally demanding and sensitive to the initialization. On the other hand, local methods are also known to be less sensitive to noise than global approaches [START_REF] Bruhn | Lucas/kanade meets horn/schunck: combining local and global optic flow methods[END_REF]. These observations suggest that the potential of local methods may still be under-exploited.

The goal of this paper is to design a new way to combine parametric models with a global variational approach through aggregation procedure, in order to both overcome the above mentioned limitations of global methods and exploit the potential of parametric estimation.

Our contributions

We propose a novel aggregation approach for optical flow estimation based on motion reconstruction from spatially varying candidates computed with parametric models.

Our method is composed of a first step estimating a collection of parametric motion models generating local motion candidates, followed by an aggregation step combining the candidates to create a global dense motion field. The main contribution of the present work is in the aggregation step. We formulate the problem as a motion reconstruction step selecting the best candidate while ensuring global smoothness of the motion field. This approach differs from other motion estimation techniques, since it decouples motion estimation and motion reconstruction. The main interest is that the reconstructed motion field is not involved in a brightness conservation constraint, and is thus not affected by its limitations. In particular, our method is able to handle large displacements without coarse-tofine schemes, it provides a valid data constraint in occluded regions, and it is more robust to noise in input images than standard variational approaches.

To achieve this, we provide motion candidates in the first step of our method that also handle large displacements, occlusions and noise in input images, by following the idea of our previous work [START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF]. We rely on the computation of parametric motion models over a set of overlapping size-variable square patches, that allows us to deal with various configurations of piecewise affine motions. An exemplar-based candidates extension strategy finds relevant motion candidates in occluded regions.

We provide an extensive experimental evaluation of our aggregation framework insisting on the versatility of its performance. We demonstrate that it outperforms the standard variational approach in case of large displacements, large occlusions and noise in input images, but also in more common situations as they can be found in the classical Middlebury benchmark. We also compare our variational aggregation with the aggregation based on discrete optimization we described in [START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF], removing any other specific features of [START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF] for fair comparison. We show that the method presented in this paper is faster and more robust to suboptimal candidate sets, while being competitive in terms of quantitative error. A first shorter version of this work was described in [START_REF] Fortun | Sparse aggregation framework for optical flow estimation[END_REF]. Compared to [START_REF] Fortun | Sparse aggregation framework for optical flow estimation[END_REF], we have integrated an occlusion handling module in the candidates estimation stage, we have modified the aggregation model to enforce the selection of a single candidate, we have improved the optimization step of our method, and we have extended the experimental validation of the method.

Related work

In this Section, we offer a brief overview of the main open issues in optical flow estimation. A recent comprehensive survey is available in [START_REF] Fortun | Optical flow modeling and computation: A survey[END_REF].

Numerous modifications of the Horn & Schunck model [START_REF] Horn | Determining optical flow[END_REF], starting with [START_REF] Black | A framework for the robust estimation of optical flow[END_REF][START_REF] Heitz | Multimodal estimation of discontinuous optical flow using markov random fields[END_REF][START_REF] Nagel | An investigation of smoothness constraints for the estimation of displacement vector fields from image sequences[END_REF], have been proposed over years, specifically to cope with large displacements and preservation of motion discontinuities [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Mémin | Dense estimation and object-based segmentation of the optical flow with robust techniques[END_REF][START_REF] Sun | A quantitative analysis of current practices in optical flow estimation and the principles behind them[END_REF][START_REF] Wedel | An improved algorithm for tv-l 1 optical flow[END_REF][START_REF] Xu | Motion detail preserving optical flow estimation[END_REF][START_REF] Zimmer | Optic flow in harmony[END_REF]. The most common response to face these two issues has been to design a multi-resolution and incremental coarse-to-fine framework along with piecewise smoothing or robust estimation. As for the data term of the global energy function other image features have been introduced like image gradient [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF], texture component [START_REF] Wedel | An improved algorithm for tv-l 1 optical flow[END_REF]. Besides, invariance properties have been sought to overcome limitations of the classical intensity constancy assumption by using Normalized Cross Correlation (NCC) [START_REF] Werlberger | Motion estimation with non-local total variation regularization[END_REF], Census transform [START_REF] Hafner | Why is the census transform good for robust optic flow computation?[END_REF], or LDP (Local Directional Pattern) descriptor [START_REF] Mohamed | Illumination-robust optical flow approach using local directional pattern[END_REF]. However, optimization complexity increases with the sophistication of the modeling.

Local and global methods may involve parametric motion models [START_REF] Black | The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields[END_REF][START_REF] Corpetti | Stochastic uncertainty models for the luminance consistency assumption[END_REF][START_REF] Cremers | Motion competition: A variational approach to piecewise parametric motion segmentation[END_REF][START_REF] Fortun | Aggregation of patch-based estimations for illumination-invariant optical flow in live cell imaging[END_REF][START_REF] Hornacek | Highly overparameterized optical flow using patchmatch belief propagation[END_REF][START_REF] Leordeanu | Locally affine sparse-to-dense matching for motion and occlusion estimation[END_REF][START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF][START_REF] Mémin | Dense estimation and object-based segmentation of the optical flow with robust techniques[END_REF][START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF][START_REF] Sun | Layered segmentation and optical flow estimation over time[END_REF][START_REF] Yang | Dense, accurate optical flow estimation with piecewise parametric model[END_REF]. The most frequently adopted ones are polynomial motion models such as translation, affine, quadratic, but other models can be investigated as well [START_REF] Hornacek | Highly overparameterized optical flow using patchmatch belief propagation[END_REF]. When attached to local optimization, the parametric motion models are estimated on local regions usually defined as square patches centered on each pixel [START_REF] Black | The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields[END_REF][START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF], possibly with an adaptation of the patch size [START_REF] Maurizot | Determination of singular points in 2D deformable flow fields[END_REF][START_REF] Senst | Robust local optical flow for feature tracking[END_REF], or its position [START_REF] Jodoin | Optical-flow based on an edge-avoidance procedure[END_REF]. This local optimization setting is easy to implement with a low computational cost, but it is clearly outperformed by sophisticated extensions of [START_REF] Horn | Determining optical flow[END_REF] in recent optical flow benchmarks [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF][START_REF] Butler | A naturalistic open source movie for optical flow evaluation[END_REF]. The motion candidates produced by our method are composed of affine motion vectors estimated in square patches without any motion segmentation. Our method implicitly selects the best patch size and position when selecting motion candidates to recover the global flow field in the second step.

When dealing with large displacements, using discrete optimization is a way to avoid resorting to coarseto-fine schemes [START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF][START_REF] Menze | Discrete optimization for optical flow[END_REF][START_REF] Yang | Dense, accurate optical flow estimation with piecewise parametric model[END_REF]. Another common approach is to somehow integrate feature correspondences in dense motion estimation. A first category of variational methods [START_REF] Braux-Zin | A general dense image matching framework combining direct and featurebased costs[END_REF][START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF][START_REF] Weinzaepfel | Deepflow: Large displacement optical flow with deep matching[END_REF] includes an additional term in the global energy. This term makes the estimated flow be close to pre-computed correspondences. However, this approach may be sensitive to matching errors by giving a fixed weight to the correspondence fitting. To overcome this problem, recent works [START_REF] Braux-Zin | A general dense image matching framework combining direct and featurebased costs[END_REF][START_REF] Timofte | Sparse flow: Sparse matching for small to large displacement optical flow[END_REF][START_REF] Weinzaepfel | Deepflow: Large displacement optical flow with deep matching[END_REF] have deliberately focused on improving the matching step. Another class of methods use correspondences to provide a coarse initialization for subsequent refinement [START_REF] Bailer | Flow fields: Dense correspondence fields for highly accurate large displacement optical flow estimation[END_REF][START_REF] Bao | Fast edge-preserving patchmatch for large displacement optical flow[END_REF][START_REF] Chen | Large displacement optical flow from nearest neighbor fields[END_REF][START_REF] Mozerov | Constrained optical flow estimation as a matching problem[END_REF][START_REF] Xu | Motion detail preserving optical flow estimation[END_REF]. In that vein, recovering a dense flow from initial sparse correspondences is also currently investigated [START_REF] Revaud | Epicflow: Edge-preserving interpolation of correspondences for optical flow[END_REF][START_REF] Timofte | Sparse flow: Sparse matching for small to large displacement optical flow[END_REF]. In [START_REF] Steinbrucker | Advanced data terms for variational optic flow estimation[END_REF], the variational refinement process is iterative and interpreted as the minimization of the orginal nonlinearized energy. The main motivation to incorporate feature matching in global optical flow methods is to alleviate the drawbacks of the coarse-to-fine scheme imposed by the classical variational optimization, in particular the loss of large displacements of small objects. Our patch correspondence substep is only involved in the motion candidates generation process and it does not drive the global optimization subsequent step.

Occlusion is a key issue in motion estimation [START_REF] Stein | Occlusion boundaries from motion: Low-level detection and mid-level reasoning[END_REF], especially in case of large displacements, since no motion measurements are available in occluded areas. By definition, a point of the current image which is occluded in the consecutive image has no corresponding point. One has to distinguish between occlusion detection, and occlusion filling with motion vectors. The two tasks can be addressed jointly within an alternate optimization strategy [START_REF] Ayvaci | Sparse occlusion detection with optical flow[END_REF][START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF][START_REF] Ince | Occlusion-aware optical flow estimation[END_REF][START_REF] Papadakis | High-dimension multilabel problems: Convex or nonconvex relaxation?[END_REF][START_REF] Sun | Symmetric stereo matching for occlusion handling[END_REF][START_REF] Sun | Local layering for joint motion estimation and occlusion detection[END_REF]. Filling occluded regions with velocity vectors given the occlusion map (or in other words, motion inpainting in occluded regions) can be related to the image inpainting problem. Image inpainting methods can be coarsely divided into diffusion-based methods [START_REF] Bertalmio | Image inpainting[END_REF][START_REF] Chan | Euler's elastica and curvature-based inpainting[END_REF] and exemplar-based methods [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF][START_REF] Komodakis | Image completion using efficient belief propagation via priority scheduling and dynamic pruning[END_REF]. Exemplar-based image inpainting fills missing parts by copying pixels of the observed image. In motion estimation, occlusion filling is usually solved by diffusion-based (or geometry-oriented) schemes, propagating motion from non-occluded regions to occluded regions using partial derivative equation (PDE) resolution [START_REF] Ayvaci | Sparse occlusion detection with optical flow[END_REF][START_REF] Berkels | Reconstructing optical flow fields by motion inpainting[END_REF][START_REF] Ince | Occlusion-aware optical flow estimation[END_REF][START_REF] Kondermann | A statistical confidence measure for optical flows[END_REF][START_REF] Papadakis | High-dimension multilabel problems: Convex or nonconvex relaxation?[END_REF][START_REF] Xu | Motion detail preserving optical flow estimation[END_REF]. In contrast, we adopt an exemplarbased strategy for candidates computation in occluded regions.

Our method share similarities with dictionary-based methods, looking for sparse combination of candidate motion vectors. Sparse representations of motion fields have recently been exploited for the design of regularization terms [START_REF] Chen | Decomposing and regularizing sparse/non-sparse components for motion field estimation[END_REF][START_REF] Dong | Nonlocal sparse and low-rank regularization for optical flow estimation[END_REF][START_REF] Jia | Optical flow estimation using learned sparse model[END_REF][START_REF] Shen | Sparsity model for robust optical flow estimation at motion discontinuities[END_REF], replacing classical spatial regularization by a proximity constraint to a sparse combination of learned patch flow fields. These strategies only act on the regularization term and are thus affected by all the above mentioned issues of global methods. Estimating directly the motion field as a linear combination of learned motion models in patches has been investigated in [START_REF] Black | Recognizing facial expressions in image sequences using local parameterized models of image motion[END_REF][START_REF] Fleet | Design and use of linear models for image motion analysis[END_REF][START_REF] Nieuwenhuis | Complex motion models for simple optical flow estimation[END_REF] with PCA decomposition on various types of training sets. However, this approach tends to produce blurry results, and has been combined with a layered approach in [START_REF] Wulff | Efficient sparse-to-dense optical flow estimation using a learned basis and layers[END_REF] to yield sharper results. One limitation is that the coefficients are estimated with a standard data term based on brightness constancy assumption. Finally, in [START_REF] Alba | Optical flow estimation with prior models obtained from phase correlation[END_REF], a pixel wise dictionary is learned online with phase correlation and a constraint on the entropy of the weights is imposed. However, the estimation only provides pixelic accuracy, without global regularization of the motion field, which causes large errors.

Robustness to noise in input images has only received little attention in the optical flow literature. In the local parametric estimation framework, explicit modeling of noise has led to dedicated methods [START_REF] Fermüller | The statistics of optical flow[END_REF][START_REF] Simoncelli | Probability distributions of optical flow[END_REF]. Experimental comparisons between local and global approaches [START_REF] Barron | Evaluation of optical flow[END_REF][START_REF] Galvin | Recovering motion fields: An evaluation of eight optical flow algorithms[END_REF] have demonstrated the highest sensitivity to noise of global approaches. Improving robustness to noise of global variational methods has been achieved in [START_REF] Bruhn | Lucas/kanade meets horn/schunck: combining local and global optic flow methods[END_REF] by integrating the local parametric assumption in the data term. However, this improvement comes at the cost of a loss of accuracy in the absence of noise.

Finally, we mention that a similar combination of candidates has been explored in the domains of image colorization [START_REF] Bugeau | Variational exemplarbased image colorization[END_REF][START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF] and image completion [START_REF] Arias | A variational framework for exemplar-based image inpainting[END_REF].

Paper organization

The rest of the paper is organized as follows. In Section 2, we present the parametric estimation of motion candidates. In Section 3, we propose an aggregation method in a variational setting. In Section 4, we demonstrate the performance of our estimation algorithms on sequences of the Middlebury and MPI Sintel datasets and other real images. Section 5 contains concluding remarks and future work.

Notations The Euclidean norm ( 2 norm) of a vector z = (z 1 , • • • , z d ) T ∈ R d is given by z 2 = ( d i=1 z 2 i ) 1/2 and the 1 norm of z by z 1 = d i=1 |z i |. The supre- mum norm of z is z ∞ = sup 1≤i≤d |z i |.
We denote two consecutive 2D image frames as u, v : Ω → R, with Ω ∈ R 2 denoting the image domain. We denote x, x or y one pixel of the image grid Ω and card(Ω) is the number of pixels.

We denote p u (x p , h)

:= u(x p + τ ), τ ∈ { h-1 2 , • • • , h+1 2 } 2 a patch of u centered at location x p ∈ Ω. The square window 1 U p (x p , h) = {x ∈ Ω : x -x p ∞ ≤ h}
is the patch support centered at pixel x p and the number of pixels falling in U p (x p , h) ⊂ Ω is h × h. We define P u := {p u (x p , h) : x p ∈ Ω, h ∈ H} as the set of all overlapping patches and H = {h 1 , . . . , h M } is a finite set of M prescribed patch sizes h m ∈ Z + .

We denote w(x) = (v 1 (x), v 2 (x)) the motion vector at pixel x of the motion field w.

The occlusion map o : Ω → {0, 1} is defined such

that o(x) = 1[x is occluded] where 1[•] is the indicator function. The set of occluded pixels is denoted O = {x ∈ Ω : o(x) = 1}.
Additional notations will be introduced in the text.

Local motion candidates and occlusion cues

We describe in this section the first step of our aggregation method. It follows the approach of [START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF] but its presentation is partly revisited. It exploits local information to supply motion candidates at each pixel. A set of motion vector candidates is generated at every pixel by a combination of patch correspondences and local parametric motion model estimations. A specific treatment is applied to occluded regions by exemplar-based extension of the motion candidates set. Our approach can be viewed as a new way to address the problem of choosing the local neighborhood for parametric estimation.

Local parametric motion candidates

The local supports for motion candidates computation are overlapping square patches of different sizes. To capture different motion scales, the patch sizes must cover a range of values. Due to the overlap and the number of patch sizes, one given pixel x ∈ Ω belongs to several patches. The candidate motion vectors at each pixel x are computed independently in each patch in two substeps described below: patch correspondences and affine motion refinement.

Patch correspondences for large displacements

We assign to each patch

p u (x p , h) in u the set {p v (y 1 , h), • • • , p v (y K , h)} of the K patches p v (•, h) in v most similar to p u (x p , h).
Hence, for each established pair of corresponding patches, we get the translation vector

t k (x p , h) ∈ Z 2 , shifting p u (x p , h) onto p v (x p + t k (x p , h) , h), k ∈ {1, • • • , K}.
Let us put forward that we do not aim at keeping at this stage the best correspondence only but at selecting K relevant correspondences to subsequently constitute motion candidates (K is assumed to be constant for all patches). The matching step is generic and could be achieved with any arbitrary feature matching algorithm (e.g., Patch-Match algorithm [START_REF] Barnes | Patchmatch: a randomized correspondence algorithm for structural image editing[END_REF]).

Affine motion refinement

The displacements estimated by patch correspondences are integer-pixel translational approximations. To attain subpixel accuracy and to allow for more complex motion, we refine the first sub-step with the estimation of a local affine motion model in every pre-registered patch pair. Denoting U p (x p , h) the pixel support of p u (x p , h), we estimate the affine motion model between two corresponding patches p u (x p , h) and p v (x p +t(x p , h), h) at pixel x = (x, y) T ∈ U p (x p , h) defined as:

δw p (x, θ p (x p , h)) = θ 1 + θ 2 x + θ 3 y θ 4 + θ 5 x + θ 6 y . ( 2
)
Assuming brightness constancy, an estimation of the parameter vector

θ p (x p , h) = (θ 1 , • • • , θ 6 ) T is the min- imizer of Up(xp,h) ϕ(v(x + δw p (x, θ p (x p , h)) + t(x p , h)) -u(x))dx (3)
where the penalty function ϕ(•) is a robust function of the family of M-estimators (e.g., Tukey's function).

Final set of motion vector candidates

The above described two-step estimation is repeated for every patch p u (x p , h) and generates a set of candidate motion vectors C(x) at each pixel x ∈ Ω. In this paper, we consider sets of regularly spaced patches, defined by a set of sizes H and an overlap ratio r ∈ [0, 1] defining the proportion of area shared by two neighbour patches of the same size. Denoting Ω p ⊂ Ω the set of center pixels of the previously defined patches, the candidates are defined as follows:

C(x) = {t k (x p , h) + δw p (x, θ p (x p , h)) : h ∈ H, (4)
x p ∈ Ω p : x p -x ∞ ≤ h, k ∈ {1, • • • , K}}.
The interest of the local set of motion candidates is first that the correspondence sub-step efficiently copes with large displacements. Specifically, it allows us to correctly deal with small structures undergoing large displacements. Second, by considering a variety of patches, we override the predefined choice of the local neighborhood. The implicit selection of the proper patch via its corresponding motion candidate is transferred to the aggregation stage. Third, introducing patches of several sizes enables to tackle motion of different scales.

Exemplar-based candidates extension in occluded regions

The computation of motion candidates described in Section 2.1 does not distinguish occluded and non-occluded pixels. However, in large occluded regions where the patches contain mostly occlusions, there is no chance to estimate relevant candidates with this local approach. Therefore, the occluded pixels require a dedicated process to compute additional motions candidates. This computation could nevertheless be considered as optional for small displacements. Indeed, considering large patch sizes enables to cope with small occlusion areas and to generate relevant candidates at motion discontinuities or at occluded positions.

When the occluded regions are known or given by an occlusion detector [START_REF] Heitz | Multimodal estimation of discontinuous optical flow using markov random fields[END_REF][START_REF] Humayun | Learning to find occlusion regions[END_REF][START_REF] Xu | Motion detail preserving optical flow estimation[END_REF], occlusion filling with motion vectors is conceptually closely related to image inpainting, since it recovers motion in regions where motion is by definition not observable. In order to deal with large occlusions produced by large displacements, we follow the inpainting analogy. In the first step of our aggregation method, the motion candidates set is thus augmented by "copy-paste" operations as described below.

We rely on the assumption that the motion at an occluded pixel x o ∈ O is similar to the motion at a close non-occluded pixel in Ω\O belonging to the same object or the same background part. The idea is to assign the set C(x) of the most similar pixel x * ∈ Ω\O to the occluded pixel x o . We limit the search for x * in a band ∂O along the occlusion boundaries. Figure 1 illustrates the matching process and the definition of O and ∂O in a simple synthetic example. Searching for the most similar pixel denoted x * ∈ Ω\O to x o is actually easier for motion inpainting than for image inpainting. Indeed, the information supplied by image u is available even in O. Thus, as x o is expected to belong to the same object as x * , we use patch similarity to find the best match in u.

An extended candidate set C + (x o ) is created for occluded pixels x o by adding to the initial set C(x o ) the motion candidates of their matched pixel x * :

C F (x o ) = C(x o ) ∪ C(x * ), ∀x o ∈ O.
(5) By convention, ∀x ∈ Ω\O, C F (x) := C(x).

A particular class of occluded (or disappearing) regions occurs at image borders in the case of large camera motion. We cope with this issue by estimating the dominant image motion due to camera motion using the Motion2D software applied to the whole image [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF], which provides additional motion candidates.

Best candidate flow

To validate our method for computing motion candidates, we have exploited sequences from MPI Sintel and Middlebury datasets [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF][START_REF] Butler | A naturalistic open source movie for optical flow evaluation[END_REF] provided with ground truth. We introduce the so-called Best Candidate Flow (BCF) by selecting at each pixel x the candidate motion vector of C F (x) closest to the ground truth vector at x. We distinguish between the BCF determined with the candidates extension described in the preceding section (or full BCF) and the BCF without it (or BCF w/o extension).

In Table 1, we report the objective evaluation given by the Endpoint Error (EPE) scores for the full BCF and BCF without candidate extensions, on the training sequences of the datasets MPI Sintel and Middlebury. Overall, the full BCF is very close to the ground truth motion field demonstrating the performance of the local parametric motion computation. We also compare these results with those of motion fields supplied by [START_REF] Weinzaepfel | Deepflow: Large displacement optical flow with deep matching[END_REF][START_REF] Xu | Motion detail preserving optical flow estimation[END_REF], as obtained with publicly available code. Clearly, full BCF outperforms these state-of-the-art methods in the two benchmarks. Accuracy is especially improved with full BCF for the MPI Sintel sequences where large displacements and wide occluded regions are present. It demonstrates that the combination of local affine estimations in square patches with patch correspondences as described in Section 2.1, is quite relevant to recover very accurate motion vectors.

Variational motion reconstruction framework

We have now to recover the global dense motion field by aggregating motion candidates available at each pixel.

We define an aggregation strategy in a variational setting, which consists in minimizing an energy of the form

E(w) = Ω ρ(w(x), C F (x)) + λ 1 φ(∇w(x))dx, (6)
where ρ(w(x), C F (x)) is a fidelity term and the second term imposes smoothness of w, balanced by the parameter λ 1 . In the following, we consider a total variation (TV) regularization: φ(∇w(x)) = ∇w(x) 1 . Unlike usual approaches for optical flow, the image intensities are not used as input of the data potential ρ(w(x), C F (x)), but are replaced by the motion candidate set C F (x). We detail in this section the modeling and optimization issues related to this reconstruction term, and the solution we adopted.

Preliminary observations

Candidates distribution

As a first investigation, we explore the information carried by the distribution of the candidates of each pixel. This analysis is motivated by the analogy with practices in other domains like image denoising or completion, where distribution of candidate image patches is exploited [START_REF] He | Image completion approaches using the statistics of similar patches[END_REF][START_REF] Salmon | Patch reprojections for nonlocal methods[END_REF]. We provide in Figure 2 six representative examples of the main forms of candidate distributions that occur in practice, and their relations with the ground truth motion vectors and the original image data. The motion vector candidates are represented by blue circles, the ground truth is the red rectangle, and the estimated motion vector is the green triangle (the full estimated motion field is given in Figure 7). In the background of the distributions, we display the value of the displaced frame difference (DFD) penalized by the 1 norm, which can be seen as a data fitting term: for the distribution at a given pixel x i (one of the six pixels in Figure 2), the value displayed in background at coordinates w

= (v 1 , v 2 ) is |v(x i + w) -u(x i )|.
The question is then to identify characteristic patterns that can allow us to identify the ground truth, given the candidates and the DFD values. Firstly, we observe that the form of the candidate distribution is highly variable. In some situations, e.g. in Figures 2(a) and 2(d), a unique mode can be clearly identified and gives a good estimate of the ground truth. However, other examples show that the main modes do not always correspond to the ground truth motion vector, and that the distribution can have highly multimodal and complex shape. In general, the form of the distribution cannot be accurately predicted from the input data. Thus, it turns out that the estimation cannot be only driven by local empirical distributions. Options like dense linear combination of candidates, fitting of a statistical distribution or clustering approaches are then not recommended.

Secondly, the relation between the DFD and the true motion vector also does not follow a general rule. It can constitute a relevant information to disambiguate complex distributions, as in Figures 2(c) and 2(f), where the ground-truth motion vector falls in regions with low values of the DFD (dark values in the background of the distribution). However, following the lowest DFD can sometimes be misleading, as it is in the case in the other figures, and it cannot be used as a unique estimation criterion.

Thirdly, cases where the true motion vector cannot be retrieved from the distribution or from low values of the DFD also occur in practice, and are illustrated by Figures 2(b) and 2(e). To handle this situation, a third information must be introduced, which can take the form of an a priori smoothness assumption.

To summarize, we have identified three sources of information to guide the design of the aggregation model: the candidate distribution, a data fitting constraint, and a smoothness assumption. These constraints are complementary and only valid locally. They should be incorporated jointly in the aggregation model in a spatially adaptive way.

Minimum distance

In addition to the qualitative analysis of the modeling aspects of the aggregation given in the previous section, another requirement that we derived from the analysis of the BCF in Section 2.3 is the selection of a single candidate at each pixel. To achieve this goal, we could define ρ(w(x), C F (x)) as the distance to the closest element of C F (x):

ρ min (w(x), C F (x)) = min i∈{1,••• ,M (x)} w(x) -w i (x) p p , (7)
where w i (x) is a motion candidate, M (x) is the number of candidates at pixel x, and p ∈ {1, 2}. The min function naturally selects one candidate used for distance measure. The proximal operator of ρ(w(x), C F (x)) can be computed exactly and the resulting energy can then be minimized in a proximal splitting framework [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF].

However, the problem of potential [START_REF] Barnes | Patchmatch: a randomized correspondence algorithm for structural image editing[END_REF] lies in its high non-convexity, leading inevitably to local minima. In practice, we experimentally observe that the algorithm converges but stays trapped in a local minimum which is very dependent on the initialization. Thus, we have to design a model that enforces the selection of a single candidate while relaxing the non-convexity of the min function [START_REF] Barnes | Patchmatch: a randomized correspondence algorithm for structural image editing[END_REF] to facilitate minimization.

Aggregation model

To this end, we introduce an additional variable α(x) = {α i (x)} i=1,..,M (x) , weighting the contribution of each candidate. The fidelity term is then expressed as

ρ(w(x), C F (x), α(x)) = M (x) i=1 α i (x) w(x) -w i (x) p p . (8)
To ensure that only one candidate is selected, the weight vector α(x) should be constrained to have binary values with only one non-zero element. To achieve this goal, we follow [START_REF] Pierre | Hue constrained image colorization in the rgb space[END_REF][START_REF] Pierre | Luminance-chrominance model for image colorization[END_REF] and point the following property: if the problem arg min i w(x) -w i (x) p p has a unique solution î, then the solution of the problem min

α(x) M (x) i=1 α i (x) w(x) -w i (x) p p , (9) s.t.
M (x) [START_REF] Barnes | Patchmatch: a randomized correspondence algorithm for structural image editing[END_REF], and is attained for α î(x) = 1 and α j (x) = 0, ∀j = î. The case where several coefficients are non-zero can only occur if the solution of arg min i w(x) -w i (x) p p is a non-singleton set S. In that case, the non-zero coefficients are {α i } i∈S and can take any configuration satisfying the constraints of [START_REF] Berkels | Reconstructing optical flow fields by motion inpainting[END_REF]. We observed that this situation rarely occurs in practice. The formulation ( 9) is convex w.r.t. to w and thus offers an algorithmically tractable alternative to the min function, while reproducing its behavior.

i=1 α i (x) = 1 ∀i ∈ {1, . . . , M (x)}, α i (x) ≥ 0, is ρ min (w(x), C F (x)) defined in
The fidelity term (8) relies only on the candidate distribution to guide the selection of a candidate. As mentioned in Section 3.1.1, purely distribution-driven estimation is insufficient to handle certain situations and should be complemented with a data-driven constraint. We exploit pre-computed confidence measures β i (x) associated to each candidate w i (x). The fidelity term is then enriched by defining ρ(w(x), C F (x), α(x)) = ( 10)

M (x) i=1 α i (x) w(x) -w i (x) p p + λ 2 β i (x) .
where λ 2 > 0 is a balance parameter. The confidence measure reflects a feature constancy assumption, e.g. based on the DFD analysed in Figure 2. Low values of β i (x) correspond to high confidence and promote high value of α i (x), such that the similarity to a distribution mode imposed by w(x) -w i (x) p p is balanced with a data fitting constraint imposed by the confidence term. Apart from [START_REF] Kondermann | A statistical confidence measure for optical flows[END_REF][START_REF] Kybic | Bootstrap optical flow confidence and uncertainty measure[END_REF], existing confidence measures are dedicated to specific motion estimation methods. For a variational approach, [START_REF] Bruhn | A confidence measure for variational optic flow methods[END_REF] uses the inverse of the global energy. For local approaches like [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF], eigenvalues of the structure tensor are usually exploited [START_REF] Mota | Analytic solutions for multiple motions[END_REF]. For parametric estimations in general, the variance of the estimate is also a possible confidence measure. To keep the generality and simplicity of our method, we consider the following simple weights based on a filtering of the DFD:

β i (x) = 1 Z Ω g(x -y)|v(y + w i (x)) -u(y)|dy, ( 11
)
where g is a convolution kernel and Z = M (x) j=1

Ω g(xy)|v(y + w j (x)) -u(y)|dy.

The analysis of Section 3.1.1 also revealed the necessity to introduce a smoothness assumption on the motion field. We complete the model with a standard Total Variation regularization to come up with the final optimisation problem min

w,α Ω N i=1 α i (x) w(x) -w i (x) p p + λ 2 β i (x) +λ 1 ∇w(x) 1 dx , (12) s.t.
M (x)

i=1 α i (x) = 1 ∀i ∈ {1, . . . , M (x)}, α i (x) ≥ 0.
This model fulfills the modeling criteria identified in Section 3.1. Minimization w.r.t. α enforces the selection of a single candidate at each pixel. The three terms of (12) combine similarity to the distribution, data-driven constraint and smoothness assumption. A key advantage of this formulation is that, differently from usual approaches based on non-linear feature conservation assumption, the optimization problem ( 12) can be solved without any linearization w.r.t w. As a result, it does not impose coarse-to-fine optimization strategies with successive linearizations at each level. Moreover, if good motion candidates have been found at occluded pixels (see Section 2.2), this data term provides a valid measure even at occlusions. It is worth noting that p = 1 enables more deviations from the candidates in case of lack of good candidates or locally wrong candidate selection.

In [START_REF] Fortun | Sparse aggregation framework for optical flow estimation[END_REF], we proposed a related model in a sparse representation framework, where the number of selected candidates was controlled by a sparsity constraint on α. The confidence measures were associated to the sparsity constraint with a weighted 1 penalization function. The drawback of this approach comes from the coupling of the sparsity constraint and the confidence measures: to ensure the selection of a single candidate, the parameter weighting the sparsity constraint has to be very high, which also gives large weight to the confidence measures. As a result, the weighted 1 term becomes predominant and the estimation is mainly driven by the confidence measures. In the model [START_REF] Black | The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields[END_REF], the selection of a candidate is decoupled from the influence of confidence measures.

Optimization

The minimization subproblems w.r.t. w and α being both convex, we resort to a block-coordinate approach alternating updates of the two variables.

Minimization w.r.t. w The minimum of (12) w.r.t. w is obtained by solving the Euler-Lagrange equations. For simplicity, we consider p = 1 in this section. The algorithm is almost equivalent for p = 2. We approximate the vectorial 1 norm by a differentiable relaxation such that z 1 ≈ ψ( z 2

2 ) = z 2 2 + 2 , with a small constant that we fix to 0.001. Under this assumption, the Euler-Lagrange equations at a given pixel x can be written:

M (x) i=1 ψ ( w(x) -w i (x) 2 2 ) α i (v j (x) -[w i (x)] j ) -λ 1 div ψ ( ∇w(x) 2 2 ) ∇v j (x)) = 0 ( 13
)
where j = {1, 2} and [•] j denotes the j th component of a vector. Using standard forward finite differences for the discretization of the gradient operator, equations ( 13) yield a non-linear system of equations, where the nonlinearity is due to the terms in ψ (•). We solve this system with the lagged nonlinearity method [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Vogel | Iterative methods for total variation denoising[END_REF]. It consits in fixing in an inner loop the nonlinear parts of ( 13), and iterating linear system solving and nonlinearity update until convergence.

Minimization w.r.t. α We solve the constrained optimization problem w.r.t. α with an Augmented Lagrangian approach [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF][START_REF] Yang | Alternating direction algorithms for \ell 1-problems in compressive sensing[END_REF]. To facilitate readability, we omit the arguments in x in this section. The positivity constraint is handled with the indicator function ι R M + defined as

ι R M + (z) = 0, if z ∈ R M + +∞, else, (14)
which leads to the following problem:

min α N i=1 α i w -w i p p + λ 2 β i + ι R M + (α), s.t. M i=1 α i = 1. ( 15
)
We reformulate ( 15) by introducing a splitting variable z associated to the indicator function:

min α,z N i=1 α i w -w i p p + λ 2 β i + ι R M + (z), s.t. z = α M i=1 α i = 1. ( 16
)
The scaled form of the Augmented Lagrangian of problem [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] writes

L(α, z, ρ 1 , ρ 2 ) = N i=1 α i w -w i p p + λ 2 β i + ι R M + (z) + µ 1 2 M i=1 α i -1 + ρ 1 µ 1 2 2 + µ 2 2 -α + z + ρ 2 µ 2 2 2 , (17)
where ρ 1 ∈ R, ρ 2 ∈ R M are Lagrange multipliers and µ 1 , µ 2 are positive penalty parameters. We use the alternated direction method of multipliers (ADMM), which separates optimization subproblems w.r.t. each variable to converge to the solution of the original problem [START_REF] Bouguet | Pyramidal implementation of the affine lucas-kanade feature tracker description of the algorithm[END_REF]. Each iteration k is composed of the following steps:

α k+1 = arg min α L(α, z k , ρ k 1 , ρ k 2 ), ( 18
)
z k+1 = arg min z L(α k+1 , z, ρ k 1 , ρ k 2 ), ( 19
)
ρ k+1 1 = ρ k 1 + µ 1 N i=1 α k+1 i -1 , (20)
ρ k+1 2 = ρ k 2 + µ 2 (-α k+1 + z k+1 ). ( 21
)
Minimization problems ( 18) and ( 19) have analytical and efficiently computable solutions. The solution of ( 18) is given by

α k+1 = µ 2 I + µ 1 1 M 1 M -1 (22
)

µ 1 M 1 - ρ k 1 µ 1 + µ 2 z k + ρ k 2 µ 2 -b λ2 ,
where we define the M -dimensional vector

1 M = (1, . . . , 1
z k+1 = max α k+1 - ρ k 2 µ 2 , 0 . ( 23
)
We emphasize that the positivity and normalization constraints of the original problem (12) define a convex set for which efficient projection can be computed, e.g. using [START_REF] Chen | Projection onto a simplex[END_REF] as proposed in [START_REF] Bugeau | Variational exemplarbased image colorization[END_REF]. However, we experimentally observed that the decoupling of the constraints yielding faster minimization subproblems in the augmented Lagrangian framework described above yielded similar results with a significantly lower computational time.

Experimental results

In this section, we analyze the performance of our VAFlow (Variational Aggregation for optical Flow) aggregation method. We highlight the versatility of VAFlow by dealing with various issues: large displacements, occlusions, motion discontinuities, noise in input images, and suboptimal candidates set. We also quantitatively demonstrate its superiority over local parametric methods and classical variational approaches on the Middlebury benchmark.

Experimental protocol

Evaluation metric When ground truth is available, we use the standard error metric for optical flow evaluation, which is the averaged endpoint error (EPE). It is defined as the average of euclidean distances at each pixel between the estimated motion vector and the ground truth.

Implementation details

The feature matching steps involved in the candidates computation (Sections 2.1.1 and 2.2) are implemented with the available code of the PatchMatch algorithm [START_REF] Barnes | Patchmatch: a randomized correspondence algorithm for structural image editing[END_REF] 2 . To achieve robustness to illumination changes, we consider a combination of saturation and value channels of the HSV color space, following [START_REF] Zimmer | Optic flow in harmony[END_REF]. The distance to minimize with PatchMatch is the sum of absolute differences (SAD) of patches. The distance between pixel x of image u and pixel y of image v is then defined for a patch size h by

SAD(x, y, h) = p u (x, h) -p v (y, h) 1 .
The affine motion estimation involved in the candidates computation step ( 3) is solved with the publicly available Motion2D software3 [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF], which implements a multi-resolution incremental minimization scheme based on the Iteratively Reweighted Least Squares (IRLS).

The occlusion detection required to extend the motion candidates in occluded regions (Section 2.2) is performed with a simple approach exploiting motion candidates computation. A coarse motion estimation is performed by block matching using PatchMatch with the smallest patch size (h = 15). The backward flow is then computed and a standard forward/backward consistency criterion [START_REF] Ince | Occlusion-aware optical flow estimation[END_REF][START_REF] Mozerov | Constrained optical flow estimation as a matching problem[END_REF] yields a coarse occlusion detection. More sophisticated methods could give more accurate occlusion regions and improve results.

In the optimization procedure described in Section 3.2.1, the motion field w is initialized by selecting at each pixel the motion candidate with best confidence measure [START_REF] Bigun | Multidimensional orientation estimation with applications to texture analysis and optical flow[END_REF]. The weights α are initialized by setting the weight corresponding to the best confidence measure to one and all the others to zero.

No post-processing is applied on the flow fields. The candidates sets were obtained with parameters H = {15, 35, 75}, r = 0.75, K = 2 (the typical number of candidates per pixel with these parameter values is around 100). The convolution filter g involved in the definition of the confidence measures in ( 11) is a rect function of size 5 × 5 pixels, which amounts to the SAD distance measure defined above and used in the patch matching step. The value of λ 2 is set to 15. Convergence of the ADMM optimization of α has been observed to be reached for 500 iterations. To save computational cost, we set a maximum number of 100 iterations, which has a limited impact on the final results. The penalty parameters are set to µ 1 = 10 and µ 2 = 10. The number of alternate optimizations in the global minimization of (12) was 4, for which convergence has been experimentally observed. The number of outer iterations involved in the fixed point iteration scheme of the minimization w.r.t. w is 15, and the linear systems are solved with the successive over-relaxation method [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF].

Methods exploited for comparison The candidates of

VAFlow are obtained with parametric estimations. Thus, the comparison with local parametric methods [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF][START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF] is informative about the efficiency of the aggregation step. In the following we will refer to "multiscale [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF]" as the coarse-to-fine implementation of [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF] described in [START_REF] Bouguet | Pyramidal implementation of the affine lucas-kanade feature tracker description of the algorithm[END_REF], and to "multiscale [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF]" as an extension of "multiscale [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF]" performing the robust affine estimation described in [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF] in each patch. The results of [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF] are obtained with the publicly available implementation4 , and we use the Motion2D software5 to apply the method [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF]. The method we call "block matching [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF]" mimics the candidates generation procedure of VAFlow. At each pixel, an initial block matching is performed and is followed by a parametric refinement between corresponding patches. Only the motion of the center pixel of the patch is kept.

As state-of-the-art results are achieved with global variational approaches, we also compare to the methods of [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF] and [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] providing open access softwares 6,7 , which implement TV-l 1 models with different optimization strategies. We also consider the method [START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF] and use the code made available by the authors8 . It extends [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] with an additional energy term imposing similarity to pre-computed feature matches. It aims at reducing the undesirable effects of the coarse-to-fine scheme. Current top performing methods [START_REF] Ranftl | Non-local total generalized variation of optical flow estimation[END_REF][START_REF] Revaud | Epicflow: Edge-preserving interpolation of correspondences for optical flow[END_REF][START_REF] Weinzaepfel | Deepflow: Large displacement optical flow with deep matching[END_REF][START_REF] Xu | Motion detail preserving optical flow estimation[END_REF] rely on the baseline principles of [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], on which they elaborate more sophisticated modules like efficient feature matching, or non-local regularization. In this paper, we propose a baseline version of our continuous aggregation concept, with simple block matching and TV regularization. Therefore, we compare it with methods [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] using the same basic ingredients. More sophisticated features could be integrated as well in our method to still improve results in the future. Finally, we also compare VAFlow with the discrete optimization approach we introduced in [START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF]. We remove the exemplar-based aggregation term and postprocessing of [START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF] to compare only the baseline aggregation methods. We refer to this method as Discrete Aggregation.

Results

Choice of the p norm We first point out the importance of the choice of p in the p norm promoting similarity to the selected candidate in [START_REF] Black | The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields[END_REF]. Table 2 gives the Average EPE obtained on the Middlebury dataset with ground truth for p = 2 and p = 1, and Figure 3 illustrates these results on an example. Choosing p = 1 yields robustness in the similarity constraint to the chosen candidate, such that few large differences between estimated motion vectors and the chosen candidate are allowed. This is a desirable property in case of locally wrong candidate selection. In Figure 3, the result with p = 2 contains two regions of large errors where the candidate selection was not optimal, whereas with p = 1, these outliers are properly handled. Few large errors could have a significant impact on the average EPE, as it can be seen in Table 2. In the light of these results, we will take p = 1 in the rest of the experiments. Chambolle and Pock [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] Brox and Malik [START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF] Fig. 4 Estimated motion fields with VAFlow and [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF], [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], [START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF] on the Bird sequence.

Large displacements of small objects One of the main limitations of coarse-to-fine schemes arises in case of large displacements of small objects, as illustrated on real sequences without ground truth in Figures 4,5 and 6. The results supplied by [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] are typical examples of failures due to coarse-to-fine schemes, which prevent here from satisfyingly recovering the duck head in Fig. 4, the ball in Fig. 5 and the foot in Fig. 6.

In contrast, VAFlow estimates correctly all these large displacements. In most cases, [START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF] also captures these movements, but at the same time, it generates large errors in other parts of the image. This is due to its high sensitivity to feature matching errors, which is better handled by VAFlow.

Motion details and discontinuities Motion details like the legs of the girl in Fig. 5 and the duck legs in Fig. 4 are better preserved by VAFlow compared to the three competing methods. In Fig. 6, the discontinuities of the motion field supplied by VAFlow are sharper and delineate better the leg and the foot of the football player.

Occlusion handling When the large displacements concern large parts of the image, occlusions become a promi- Chambolle and Pock [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] Brox and Malik [START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF] Fig. 5 Estimated motion fields with VAFlow and [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF], [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], [START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF] on the Backyard sequence of the Middlebury dataset.

nent issue, as illustrated in the three image pairs of Fig. 7. To demonstrate the effect of our occlusion handling, we desactivate the occlusion handling module (VAFlow w/o occlusions) in the motion candidate generation step, and compare the results with those obtained by the full VAFlow method. In each case, VAFlow w/o occlusion still captures well large displacements, but it also exhibits large errors at occluded pixels, due to the absence of good candidates. When occlusion handling is activated, the result is visually greatly improved in these regions and is very close to the ground truth. This observation is confirmed by the large decrease of the EPE in each case (also reported in Fig. 7).

Quantitative evaluation We provide a quantitative evaluation in Table 3, reporting the EPE obtained with VAFlow, local approaches [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF][START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF], and variational methods [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] for the sequences of the Middlebury dataset with ground truth. The candidates of VAFlow are computed by local methods. In particular, they are obtained with the same estimation procedure as block matching [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF]. Therefore, the large improvement offered by VAFlow on these methods is due to the efficiency of the aggregation step, which is able to select the best motion candidate rather than just keeping the motion estimate at the central point of each patch. VAFlow also outperforms the global variational approaches [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] on almost all the sequences. Chambolle and Pock [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] Brox and Malik [START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF] Fig. 6 Estimated motion fields with VAFlow and methods [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF], [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], [START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF] on the Football sequence.

In Table 4, we report results obtained on the MPI Sintel training dataset [START_REF] Butler | A naturalistic open source movie for optical flow evaluation[END_REF], characterized by the presence of sequences with very large displacements. We give the average error on the whole benchmark, and we also give average errors obtained on the seven sequences with the largest displacements. The advantage between of VAFlow over the other methods is larger than in Table 3, which confirms the ability of our aggregation strategy to handle large displacements, in particular compared to the integration of feature matching as an additional constraint in a classical variational approach [START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF]. We mention some recent methods like [4] and [START_REF] Menze | Discrete optimization for optical flow[END_REF] are able to outperform these results with an average endpoint error of respectively 2.61 and 2.25 on the whole MPI Sintel training dataset. However, as explained in Section 4.1, these methods exploit sophisticated modules that could be integrated in our framework. For instance, the contributions of [START_REF] Bailer | Flow fields: Dense correspondence fields for highly accurate large displacement optical flow estimation[END_REF] and [START_REF] Menze | Discrete optimization for optical flow[END_REF] are coarse feature matching methods that have to be refined with the variational method [START_REF] Revaud | Epicflow: Edge-preserving interpolation of correspondences for optical flow[END_REF], which integrates a sophisticated edge detection to sharpen motion discontinuities. These ingredients could be incorporated in our aggregation model to improve results and compete on state-of-the-art computer vision benchmarks. Our primary aim is to propose a general aggregation framework for motion estimation.

We also report in Table 3 and Table 4 results obtained by selecting at each pixel the candidate with the best confidence measure (the lowest value of β i (x)), which we refer to as "Best confidence flow". The results are always significantly worse than those of VAFlow. It demonstrates that the motion estimation with VAFlow is not over-guided by the confidence measures and can deviate from them to improve global accuracy of the motion field.

Values of the selection weights α The averaged final estimation of α obtained for the whole sequence temple 2 of the MPI Sintel dataset (50 frames) is illustrated in Figure 8. 98.7% of the coefficients are lower than 0.01 and are considered to have no siginificant influence on the final results. Therefore, only coefficients superior to 0.01 are displayed in Figure 8. Among coefficients greater than 0.01, 96% are greater than 0.95, which confirms that the algorithm selects only one candidate for the reconstruction most of the time. In that sense, our method finds the sparsest solution in most cases.

Robustness to noise Existing optical flow benchmarks do not integrate robustness to noise as an evaluation criterion. However, it is common to deal with noisy images when specific optical devices are used, as in microscopy or astronomy.

VAFlow performs patch-based parametric motion estimation, in the candidates generation step. The aggregation step (motion reconstruction) does not exploit any pixel wise feature conservation assumption, but only uses a patch-based confidence measure. Parametric estimations in patches [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF][START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF] are known to be more robust to noise than global variational methods. There-fore, we expect VAFlow to provide with robustness to noise while ensuring of the accuracy global variational methods in the absence of noise, as previously demonstrated in Table 3.

In Fig. 9, we plot the average EPE for Middlebury sequences with ground truth after adding Gaussian noise to the input images with different standard deviations. The results supplied by VAFlow are compared with those of [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] in Fig. 9.a and with [START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF] Table 3 Endpoints error obtained with VAFlow, the local methods [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF][START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF] and the variational methods [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF] on the Middlebury dataset with ground truth. in Fig. 9.b. The impact of noise is significantly lower on the performance of VAFlow than on those of [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. The difference is even more pronounced between VAFlow and [START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF], which is due to the high sensitivity of [START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF] to wrong feature matches, as already emphasized in previous results.

Suboptimal candidates set

The final output of VAFlow is dependent on the quality of the motion candidates. More patches should be considered to augment the variety of candidates. A crucial parameter is the overlap ratio r ∈ [0, 1], defining the amount of common area shared by two neighbor patches. When r is close to one, there are as many patches as pixels for a given patch size, and the number of candidates the is highest. However, the number of patches also increases the computation time, such that a trade-off has to be found between accuracy and complexity. The impact of the overlap ratio on these two aspects is reported in Table 5, which summarizes the evolution with r of the average EPE on the Middlebury benchmark sequences with ground truth on and the computational time on the Urban 2 sequence of the Middlebury dataset. While the computation time increases slowly when r is small, it changes much faster when r > 0.5. In the same time the error increase remains relatively limited for r > 0.5. This robustness to suboptimal candidates sets is further emphasized by the visual results of Figure 10, where we can observe that the results stay very similar when r decreases, in particular when r > 0.5. In practical scenarios where computational time matters, this robustness can allow us to make huge gains in complexity without loosing too much accuracy.

Comparision with discrete optimization [START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF] We focus now on the comparison between the variational aggregation scheme of VAFlow and the aggregation based on discrete optimization described in [START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF], that we call Discrete Aggregation. general slightly more accurate than those of VAFlow. However, the advantage of VAFlow lies in its robustness to suboptimal candidate sets and its computation time.

Conclusion

We have proposed a variational aggregation framework for optical flow estimation based on a sparse representation of the motion field. We combine in two successive steps local parametric estimation yielding motion candidates, and global aggregation supplying the global recovered flow. We formulated the aggregation step as a global energy minimization problem without coarse-tofine strategy, combining the best motion candidates at every pixel while preserving motion discontinuities. We promoted sparse solutions, that is, the selection at each pixel of a few candidates in space-variant motion vector dictionaries. We handle occlusion with an exemplarbased motion inpainting approach in the candidates computation step. We demonstrated the improvements yielded by our method over standard variational approaches in various situations of large displacements of small objects, occlusions, noise in input images and motion discontinuities. We also achieved a lower computational time and more robustness to suboptimal candidates set compared to the discrete aggregation approach introduced [START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF]. The framework is generic, and both the local and global steps could be adapted for specific purposes, especially using more sophisticated feature matching techniques.

Frame 1 Frame 2 Fig. 1

 21 Fig. 1 Illustration of the exemplar-based inpainting of motion candidates. The foreground is shifting to the right over a static background. The candidate set of occluded pixel xo ∈ O (in red) is extended by adding the candidates of its matched non-occluded pixel x * ∈ ∂O (in yellow).

Fig. 2

 2 Fig.2Visualization of the distribution of the motion candidates at several pixels in the image. The central image is the ground truth motion field of the frame 23 of the temple 2 sequence of the MPI Sintel data set. The six plots represent the motion vector candidates (blue circles), the motion vector ground truth (red rectangle) and the estimated motion vector with our method (green triangle)(the estimated motion field on the whole image is given in Figure7) at each corresponding pixel. The horizontal and vertical axes are respectively the horizontal and vertical components of the motion vectors. The value of the displaced frame difference (DFD) penalized by the 1 norm is displayed in the background of the distributions.

  ) and the vector b λ2 = (b 1 , . . . , b M ) with b i = p ww i p p +λ 2 β i . The matrix inversion can be easily achieved with the Sherman-Morrison formula. The update of z is a simple projection onto the set R M + given by

p = 2 p = 1 Average 1 Fig. 3

 113 Fig. 3 Illustration of of the impact of the choice of p for the data fidelity term.

Frame 1 Frame 2 VAFlow

 2 Brox et al.[START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF]

Fig. 8

 8 Fig. 8 Average distribution of the coefficients α on the temple 2 sequence of the MPI Sintel benchmark (50 frames). Only coefficients greater than 0.01 are displayed. 98.7% of the coefficients are below 0.01.

ambush 2 -

 2 Frame 1 ambush 2 -Frame 2 Occlusions ground truth (in white) Motion field ground truth VAFlow w/o occlusions -EPE=17.909 VAFlow -EPE=4.363 temple 3 -Frame 19 temple 3 -Frame 20 Occlusions ground truth (in white) Motion field ground truth VAFlow w/o occlusions -EPE=9.495 VAFlow -EPE=2.545 temple 2 -Frame 23 temple 2 24 -Frame 2 Occlusions ground truth (in white) Motion field ground truth VAFlow w/o occlusions -EPE=7.474 VAFlow -EPE=3.286

Fig. 7

 7 Fig. 7 Comparison of motion fields computed with VAFlow and with VAFlow without the occlusion handling module (w/o occlusions). Results are obtained on sequences of the MPI Sintel dataset [23] with large displacements. The EPE of each result is given in the captions attached to the motion fields.

Fig. 9 Fig. 10

 910 Fig.9Evolution of the EPE with the standard deviation of the added Gaussian noise in the input images. The reported EPE is the average EPE over all the sequences of the Middlebury dataset with ground truth. Fig.9acompares VAFlow with[START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] and Fig.9bcompares VAFlow with[START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF].

Figure 11

 11 compares the impact of the overlap ratio on the EPE and the computation time. While the EPE of Discrete Aggregation is lower for a large overlap ratio, the results of VAFlow are less impacted by a lower quality of the candidates, and it gives lower EPE when the overlap ratio is approximately below 0.45. In the same time, the computation time of Discrete Aggregation increases faster than the one of VAFlow with r. For r = 0.75, Discrete aggregation is almost two times slower than VAFlow.

Fig. 11

 11 Fig. 11 Comparison of the behaviour of VAFlow and Discrete Aggregation w.r.t. the overlap ratio. Fig. (a) show the evolution of the average EPE on the sequences with ground truth of the Middlebury benchmark, and Fig. (b) shows the evolution of the computational time on the Urban 2 sequence of the Middlebury benchmark.

Table 1

 1 EPE-all scores of motion fields on sequences with ground truth from MPI Sintel and Middlebury datasets

		MPI Sintel	Middlebury
	Full BCF	0.792	0.0710
	BCF w/o extension	1.851	0.0833
	DeepFlow [83]	4.691	0.386
	MDP-Flow2 [86]	4.006	0.223

Table 2

 2 Average EPE results on the Middlebury benchmark for p = 1 and p = 2.

Table 4

 4 Endpoints error obtained with VAFlow, the local methods[START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF][START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF] and the variational methods[START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF] on the MPI Sintel dataset with ground truth. The last column gives the average result on the whole data set, and results on the seven sequences with the largest displacements are also given.

	Sequences with	cave 2 market 6 temple 3 ambush 5 ambush 6 ambush 2 market 5	Average on
	large displacements								whole benchmark
	Best confidence flow	22.02	12.95	21.19	16.86	25.90	26.13	38.87	11.19
	VAFlow	7.99	4.82	8.74	6.34	7.86	10.17	11.79	3.90
	Brox et al. [18]	27.54	7.30	15.84	12.72	15.44	34.94	23.07	7.31
	Chambolle and Pock [24]	25.01	8.55	21.43	12.22	16.07	35.67	23.74	7.91
	Brox and Malik [19]	9.20	5.61	14.67	10.90	11.11	20.73	14.98	5.03

Table 5

 5 Evolution with the overlap ratio r of the Average EPE on the Middlebury dataset with ground truth and the computational time on the Urban 2 sequence of the Middlebury dataset.

	Overlap ratio	0.75	0.5	0.25	0.1
	Average EPE	0.296 0.310 0.329 0.354
	Computation time (s)	305	142	117	111

  Table 6 reports the EPE obtained on sequences of the Middlebury and MPI Sintel datasets with ground truth by VAFlow and Discrete Aggregation. Results supplied by Discrete Aggregation are in

Table 6

 6 Endpoints error obtained with VAFlow and Discrete Aggregation on the sequences with ground truth of the Middlebury and MPI Sintel datasets.

	Middleburry		Grove2 Grove3 Urban2 Urban3 Venus RubberWhale Dimetrodon Hydrangea
	VAFlow		0.161		0.630	0.374		0.395	0.298	0.134		0.090	0.194
	Discrete aggregation	0.166		0.621	0.337		0.381	0.287	0.121		0.122	0.179
	MPI Sintel		cave 2 market 6 temple 3 ambush 5 ambush 6 ambush 2 market 5
	VAFlow		7.99		4.82	8.74		6.34	7.86		10.17		11.79
	Discrete aggregation	8.228		4.547	8.314	5.50	6.251		9.456		11.958
		0.4									400
						Discrete aggregation			Discrete aggregation
						VAFlow					VAFlow
	EPE	0.3 0.32 0.34 0.36 0.38								Computational time (s)	200 250 300 350
		0.28									150
		0.1 0.26	0.2	0.3	0.4	0.5	0.6	0.7	0.8		0.1 100	0.2	0.3	0.4	0.5	0.6	0.7	0.8
				Overlap ratio								Overlap ratio

Without loss of generality, isotropic circular patches could be considered as well.

http://gfx.cs.princeton.edu/pubs/Barnes 2009 PAR/index.php

http://www.irisa.fr/vista/Motion2D/

http://www.mathworks.com/matlabcentral/fileexchange/23142iterative-pyramidal-lk-optical-flow

[START_REF] Baker | A database and evaluation methodology for optical flow[END_REF] http://www.irisa.fr/vista/Motion2D

/ 6 http://lmb.informatik.uni-freiburg.de/resources/software.

php 7 http://gpu4vision.icg.tugraz.at/index.php?

content=downloads.php 8 http://lmb.informatik.uni-freiburg.de/resources/software.php

Acknowledgements This work was realized as part of the Quaero program, funded by OSEO, French State agency for innovation. The authors acknowledge France-BioImaging infrastructure supported by the French National Research Agency (ANR-10-INBS-04-07, "Investments for the future"). They thank also the reviewers for useful comments helping improving the paper. Finally, they thank Ferreol Soulez, Martin Storath, Olivier Demetz, Simon Setzer and Joachim Weickert for inspiring discussions at different stages of this work.