Denis Fortun
email: denis.fortun@epfl.ch

Patrick Bouthemy

Charles Kervrann

Optical flow modeling and computation: a survey

Keywords: Optical flow, Motion estimation, Regularization, Parametric models, Optimization, Feature matching, Occlusions

Optical flow estimation is one of the oldest and still most active research domains in computer vision. In 35 years, many methodological concepts have been introduced and have progressively improved performances, while opening the way to new challenges. In the last decade, the growing interest in evaluation benchmarks has stimulated a great amount of work. In this paper, we propose a survey of optical flow estimation classifying the main principles elaborated during this evolution, with a particular concern given to recent developments. It is conceived as a tutorial organizing in a comprehensive framework current approaches and practices. We give insights on the motivations, interests and limitations of modeling and optimization techniques, and we highlight similarities between methods to allow for a clear understanding of their behavior.

Introduction

Motion analysis is one of the main tasks of computer vision. From an applicative viewpoint, the information brought by the dynamical behavior of observed objects or by the movement of the camera itself is a decisive element for the interpretation of observed phenomena. The motion characterizations can be extremely variable among the large number of application domains. Indeed, one can be interested in tracking objects, quantifying deformations, retrieving dominant motion, detecting abnormal behaviors, and so on. The most low-level characterization is the estimation of a dense motion field, corresponding to the displacement of each pixel, which is called optical flow. Most high-level motion analysis tasks employ optical flow as a fundamental basis upon which more semantic interpretation is built.

Optical flow estimation has given rise to a tremendous quantity of works for 35 years. If a certain continuity can be found since the seminal works of [START_REF] Horn | Determining optical flow[END_REF][START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF], a number of methodological innovations have progressively changed the field and improved performances. Evaluation benchmarks and applicative domains have followed this progress by proposing new challenges allowing methods to face more and more difficult situations in terms of motion discontinuities, large displacements, illumination changes or computational costs. Despite great advances, handling these issues in a unique method still remains an open problem. Comprehensive surveys of optical flow literature were carried out in the nineties [START_REF] Beauchemin | The computation of optical flow[END_REF][START_REF] Mitiche | Computation and analysis of image motion: A synopsis of current problems and methods[END_REF][START_REF] Stiller | Estimating motion in image sequences[END_REF]. More recently, reviewing works have focused on variational approaches [START_REF] Weickert | A survey on variational optic flow methods for small displacements[END_REF], benchmark results [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF], specific applications [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF], or tutorials restricted to a certain subset of methods [START_REF] Mitiche | Computer Vision Analysis of Image Motion by Variational Methods[END_REF][START_REF] Wedel | Optical flow estimation[END_REF].

However, covering all the main estimation approaches and including recent developments in a comprehensive classification is still lacking in the optical flow field. This survey paper is conceived as a tutorial aiming at organizing in a comprehensive framework the main approaches and practices for optical flow estimation. We will not try to tell the whole story of optical flow evolution, neither we give an exhaustive list of existing methods. We rather propose a synthetic classification of the main methodological principles of existing methods, with a particular concern given to recent developments. We will insist on the modeling aspects, practical interests and limitations of each introduced methodological element, with the aim of providing useful cues for a better understanding of optical flow. We adopt a classifying approach, decomposing the optical flow estimation problem in sub-parts. This viewpoint has the advantage of being didactic and facilitates the identification of similarities between methods in a global unique framework. However, this ambition of clarity should not hide that individual methods are often conceived as coherent approaches, and cannot be reduced to a mere assembly of isolated parts.

The remainder of this paper is organized as follows. In Section 2, we present general and basis concepts on optical flow used throughout the paper. Section 3 is devoted to a review of data constancy assumptions for optical flow and their usage for the design of data terms. In Section 4, the parametric approach, its specific issues and existing solutions are described. Section 5 presents the modelling and optimization aspects of globally regularized methods. Section 6 focuses on the occlusion issue. Section 7 is dedicated to the recent trend of integrating feature matching in the optical flow estimation process. Finally, the general conclusions of this survey are reported in Section 8.

Preliminaries

Definition of optical flow

The notion of optical flow literally refers to the displacements of intensity patterns. This definition originates from a physiological description of the visual perception of the world through image formation on the retina. In that sense, while optical flow is necessarily caused by relative motion between the observer and the objects of the observed scene, it only represents motion of intensities in the image plane, and not necessarily accounts for the actual 3D motion in the physical scene [START_REF] Verri | Motion field and optical flow: Qualitative properties[END_REF]. This definition of optical flow states that intensity of moving pixels remains constant during motion, and the main difficulty is to cope with the uncertainty due to the ill-posedness of the problem, called aperture problem (described in subsection 2.4).

However, motion of interest in computer vision is in practice the real displacements of the objects. Thus, the variable that we actually want to retrieve is the projection on the image plane of the 3D motion in the scene, usually called motion field. The problem of recovering the 3D motion itself from several synchronized cameras is also an active research domain known as scene flow estimation [START_REF] Basha | Multi-view scene flow estimation: A view centered variational approach[END_REF][START_REF] Vedula | Threedimensional scene flow[END_REF]. The new problem arising is then that intensity changes are not necessarily caused by the displacements of objects of the scene, but can also be due to other disturbing phenomena like lighting changes, reflection effects or modifications of the properties of the objects affecting their light emission or reflectance.

The terms optical flow and motion field are usually mixed up to both designate the above definition of the motion field, and we will carry on this practice in the sequel of this survey.

Applicative domains

Optical flow provides fundamental information at the basis of various computer vision systems, in a wide range of applicative domains. We cite here some of the main fields where it is currently exploited.

The increasing use of multimedia devices involving video displays accentuates the importance of several video analysis problems. Action recognition is for instance an essential task for semantic interpretation of video content [START_REF] Ali | Human action recognition in videos using kinematic features and multiple instance learning[END_REF][START_REF] Jain | Better exploiting motion for better action recognition[END_REF][START_REF] Wang | Dense trajectories and motion boundary descriptors for action recognition[END_REF]. Video compression standards like MPEG exploit motion estimation to predict intermediate frames [START_REF] Jakubowski | Block-based motion estimation algorithmsa survey[END_REF]. Optical flow can also be a determinant feature in video indexing and retrieval [START_REF] Hu | A survey on visual content-based video indexing and retrieval[END_REF][START_REF] Piriou | Recognition of dynamic video contents with global probabilistic models of visual motion[END_REF][START_REF] Su | Motion flow-based video retrieval[END_REF]. Video restoration of aged films is another field of interest, and can be improved by taking optical flow estimation into account [START_REF] Gal | Progress in the restoration of image sequences degraded by atmospheric turbulence[END_REF][START_REF] Werlberger | Optical flow guided tv-l1 video interpolation and restoration[END_REF].

In a biomedical context, dynamic properties of tissues or cellular objects are also of upmost importance. The deformation of organs [START_REF] Hata | Three-dimensional optical flow method for measurement of volumetric brain deformation from intraoperative mr images[END_REF][START_REF] Xavier | An adapted optical flow algorithm for robust quantification of cardiac wall motion from standard cine-mr examinations[END_REF] or the estimation of blood flow [START_REF] Huang | Quantification of blood flow in internal cerebral artery by optical flow method on digital subtraction angiography in comparison with time-of-flight magnetic resonance angiography[END_REF] are examples of medical applications that can require optical flow computation. Many conceptual and methodological links can also be found between medical image registration and optical flow [START_REF] Glocker | Inter and intra-modal deformable registration: Continuous deformations meet efficient optimal linear programming[END_REF][START_REF] Rueckert | Nonrigid registration using free-form deformations: application to breast MR images[END_REF][START_REF] Schnabel | A generic framework for non-rigid registration based on nonuniform multi-level free-form deformations[END_REF][START_REF] Shi | Registration using sparse free-form deformations[END_REF][START_REF] Sotiras | Deformable medical image registration: A survey[END_REF]. In microscopy, dense motion can inform about cell deformation [START_REF] Amat | Fast and robust optical flow for time-lapse microscopy using super-voxels[END_REF][START_REF] Fortun | Correlation and variational approaches for motion and diffusion estimation in fluorescence imaging[END_REF][START_REF] Kim | Nonrigid registration of 2-d and 3-d dynamic cell nuclei images for improved classification of subcellular particle motion[END_REF][START_REF] Pizarro | Towards dense motion estimation in light and electron microscopy[END_REF], motion of cellular structures [START_REF] Delpiano | Performance of optical flow techniques for motion analysis of fluorescent point signals in confocal microscopy[END_REF][START_REF] Fortun | Aggregation of patch-based estimations for illumination-invariant optical flow in live cell imaging[END_REF], or help for individual cell tracking [START_REF] Liu | Optical flow guided cell segmentation and tracking in developing tissue[END_REF]. Robots or vehicle navigation also exploit optical flow as input of control systems for automatic guidance. Autonomous car driving has received particular attention in recent years [START_REF] Geiger | Are we ready for autonomous driving? the KITTI vision benchmark suite[END_REF][START_REF] Giachetti | The use of optical flow for road navigation[END_REF][START_REF] Sun | On-road vehicle detection: A review[END_REF]. Obstacle detection and avoidance are the main tasks investigated for general robot control in real environment exploiting optical flow [START_REF] Chao | A survey of optical flow techniques for robotics navigation applications[END_REF][START_REF] Crétual | Visual servoing based on image motion[END_REF][START_REF] Enkelmann | Obstacle detection by evaluation of optical flow fields from image sequences[END_REF].

Automated video surveillance is another growing research field requiring motion analysis. In particular, facial expression and gesture recognition [START_REF] Black | Recognizing facial expressions in image sequences using local parameterized models of image motion[END_REF][START_REF] Cutler | View-based interpretation of real-time optical flow for gesture recognition[END_REF], crowd motion and pedestrian behavior analysis [START_REF] Basset | Recovery of motion patterns and dominant paths in videos of crowded scenes[END_REF][START_REF] Kiryati | Real-time abnormal motion detection in surveillance video[END_REF] often involve dense motion analysis.

Among the other fields of interest, fluid flow analysis has applications in meteorology, oceanography aerodynamics or fluid mechanics and requires to consider specific physical constraints [START_REF] Corpetti | Dense estimation of fluid flows[END_REF][START_REF] Héas | Layered estimation of atmospheric mesoscale dynamics from satellite imagery[END_REF][START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF][START_REF] Liu | Fluid flow and optical flow[END_REF]. Several typical image deformations characterized by the conservation of texture properties during motion can be described as dynamic textures, and receive a specific treatment [START_REF] Crivelli | Motion textures: Modeling, classification, and segmentation using mixed-state Markov random fields[END_REF][START_REF] Fazekas | Normal versus complete flow in dynamic texture recognition: a comparative study[END_REF][START_REF] Fazekas | Dynamic texture detection based on motion analysis[END_REF].

Evaluation of optical flow methods

As in other computer vision domains, much attention has been paid to the design of appropriate evaluation procedures for optical flow.

The visualization of motion fields provides a qualitative insight on the accuracy of the estimation. The two main visualization techniques are presented in Figure 1. The arrow visualization directly represents motion vector and offers a good intuitive Barron et al. [START_REF] Barron | Evaluation of optical flow[END_REF] Middlebury [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF] Kitti [START_REF] Geiger | Are we ready for autonomous driving? the KITTI vision benchmark suite[END_REF] MPI Sintel [START_REF] Butler | A naturalistic open source movie for optical flow evaluation[END_REF] perception of physical motion. On the counterpart, a clean display requires to under-sample the motion field to prevent overlapping of arrows. The color code visualization associates a color hue to a direction and a saturation to the magnitude of the vector. It allows for a dense visualization of the flow field and a better visual perception of subtle differences between neighbor motion vectors.

On the other hand, objective quantitative evaluation based on error metrics is necessary for an accurate comparison of method performances. When ground truth is available, two error measures are commonly used, namely the Angular Error (AE) and the Endpoint Error (EPE). The AE of an estimated motion vector w est = (u est , v est ) w.r.t. the reference vector w re f = (u re f , v re f ) is defined by the 3D angle created by the extended vectors (u est , v est , 1) and (u re f , v re f , 1) :

AE = cos -1              u est u re f + v est v re f + 1 u 2 est + v 2 est + 1 u 2 re f + v 2 re f + 1              . ( 1
)
The EPE is defined as the Euclidean distance between the two vectors:

EPE = (u est -u re f ) 2 + (v est -v re f ) 2 . (2)
These two metrics are complementary regarding motion magnitude. AE is very sensitive to small estimation errors, occurring for small displacement, whereas EPE hardly discriminates between close motion vectors. On the contrary, AE tends to under-evaluate distances between motion vectors in some cases of large motion magnitudes, while EPE strongly penalizes large estimation errors.

The design of challenging benchmarks with ground truth on which to apply these error measures has motivated a substantial amount of work, which main steps are illustrated in Figure 2. The first optical flow benchmark with ground truth established by [START_REF] Barron | Evaluation of optical flow[END_REF] dealt either with simple parametric transformations applied to real images, like translation or rotation, or with synthetic sequences for which the true motion is available by construction. The resulting motion fields were characterized by small displacements and absence of discontinuities. More recent benchmarks have been designed to address new types of problems or specific applicative issues. The Middlebury benchmark 1 [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF][START_REF] Baker | A database and evaluation methodology for optical flow[END_REF] is composed of more challenging sequences, partly made of smooth deformations similar to the sequences described in [START_REF] Barron | Evaluation of optical flow[END_REF], but also involving motion discontinuities and motion details. While some sequences are synthetic, several others were acquired in a strictly controlled environment allowing to produce ground truth for real scenes. Issues raised by Middlebury being almost solved by modern methods, the MPI Sintel benchmark 2 [START_REF] Butler | A naturalistic open source movie for optical flow evaluation[END_REF] has been recently proposed. It is extracted from a synthetic movie and opens new issues mostly related to very large displacements, occlusions, illumination 1 http://vision.middlebury.edu/flow/eval/ 2 http://sintel.is.tue.mpg.de/ changes, and effects like blur or defocus. In parallel, dedicated datasets have been designed to solve specific problems related to applicative contexts, the most successful example being the KITTI benchmark3 [START_REF] Geiger | Are we ready for autonomous driving? the KITTI vision benchmark suite[END_REF] devoted to assisted driving applications.

Estimation principles

Let us denote an image sequence by I : Ω × T → R, where Ω ⊂ R 2 is the image domain and T is the sampled time interval of the sequence. Every optical flow estimation method is based on an assumption on the relationship between the searched motion field w : Ω → R 2 at time t and the image I(•, t). The most natural and widely used assumption is that pixel intensity remains constant during displacement. The brightness constancy constraint equation (BCCE) is then defined by:

dI dt (x(t), t) = 0. ( 3
)
Other feature constancy can be chosen, each encoding specific image properties, which will be discussed in Section 3. The discrete approximation of (3) at a given pixel x ∈ Ω and time t yields:

I(x + w(x), t + 1) -I(x, t) = 0. (4)
However, constraint (4) usually generates particularly difficult optimization problems. It can be much more tractable to consider the expanded version of (3) with partial derivatives, resulting in a linear version of (4):

∂I ∂x 1 (x) u(x) + ∂I ∂x 2 (x) v(x) + ∂I ∂t (x) = 0, (5)
where w(x) = (u(x), v(x)) and x = (x 1 , x 2 ) . For the sake of clarity, we have dropped the dependency of I over t. We can also write [START_REF] Amiaz | Coarse to over-fine optical flow estimation[END_REF] as

∇I(x) • w(x) + I t (x) = 0, (6)
where

∇• = ∂• ∂x 1 , ∂• ∂x 2 T
is the spatial gradient operator, I t = ∂I/∂t is the partial temporal derivative and the dot " • " denotes the inner product.

The linearized brightness constancy constraint (5) provides only one equation to recover the two unknown components of w(x). From this single constraint, the component of the motion vector w(x) in the direction of the image gradient can be computed, but the two-dimensional problem remains underconstrained. This is known as the aperture problem, stating that motion of linear structures, as it is assumed by [START_REF] Amiaz | Coarse to over-fine optical flow estimation[END_REF], is by nature ambiguous if the neighboring context is not taken into account.

To make the problem well-posed, it is necessary to introduce an additional constraint encoding a priori information on w. The a priori will take the form of spatial coherency imposed by either local or global constraints, described in Section 4 and Section 5.

It is important to mention that while (4) holds for motion of arbitrary magnitude, the continuous motion constraint (5) restricts its validity domain to the linear region of I, which usually corresponds to small displacements or very smooth images. The linearization is nevertheless necessary for methods relying on differential computations. The standard technique to cope with large displacements is to embed the estimation in a coarse-to-fine scheme [START_REF] Black | The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields[END_REF][START_REF] Enkelmann | Investigations of multigrid algorithms for the estimation of optical flow fields in image sequences[END_REF][START_REF] Mémin | Dense estimation and object-based segmentation of the optical flow with robust techniques[END_REF]. The idea is to create a pyramid of coarse-to-fine downsampled versions of the original image. At the coarsest level, the linearity domain of the image encompasses large displacements and the estimation can be based on [START_REF] Amiaz | Coarse to over-fine optical flow estimation[END_REF]. The estimations at coarser levels serves to warp the image at subsequent finer levels, where the estimation then reduces to a search for small motion increments. The solution is iteratively refined at each level until reaching the full image resolution. It is possible to continue the process in still higher resolutions created by interpolation, as proposed in the coarse-to-overfine approach [START_REF] Amiaz | Coarse to over-fine optical flow estimation[END_REF]. The solution at each level of the multi-resolution pyramid can be interpreted as a fixed point in the direct optimization of the non-linearized constancy equation (4) [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF]. Almost all differential methods concerned with large displacements resort to the multiscale approach, possibly with some additional strategies to avoid its drawbacks.

The main undesirable effect produced by smoothing at coarse levels is the loss of small and rapidly moving objects in the final estimated flow field. If the object extent is smaller than its displacement, it is likely to be smoothed out at coarse levels and then "forgotten". Avoiding this drawback has been a very active topic in recent years. It has been accomplished mainly by integrating feature matching in the estimation process, as will be discussed in Section 7, or by resorting to discrete optimization methods, as detailed in Section 5.2.2.

Data cost

The deviations from the data constancy assumption (3) are penalized at every pixel x ∈ Ω via a data potential ρ data (x, I 1 , I 2 , w), defined in the case of brightness constancy (4) as

ρ data (x, I 1 , I 2 , w) = φ(I 2 (x + w(x)) -I 1 (x)), (7)
where I 1 = I(•, t) and I 2 = I(•, t + 1) denote two successive frames, and φ(•) is the penalty function.

The brightness constancy assumption is in practice an imperfect photometric expression of the real physical motion in the scene. A common counter-example consists in moving the light source of an immobile scene, producing brightness variations without motion of any objects. In general, while it is possible to create synthetic sequences for which the constraint strictly holds [START_REF] Butler | A naturalistic open source movie for optical flow evaluation[END_REF], it is often violated in practice in case of changes in the illumination sources of the scene, shadows, noise in the acquisition process, specular reflections or large and complex deformation.

Choosing a quadratic penalty function φ(z) = z 2 , as in early works [START_REF] Horn | Determining optical flow[END_REF][START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF], makes optimization much easier, but assumes that the residual of the brightness constancy constraint equation (3) is normally distributed and thus gives a strong influence to large localized violations mentioned above. It is then common to resort to robust statistics [START_REF] Huber | Robust statistics[END_REF] to reduce the impact of local errors considered as outliers. Common robust alternatives to quadratic penalty are the L 1 norm [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF], the Tukey function [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF], the Lorentzian norm [START_REF] Black | The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields[END_REF] or the Leclerc's function [START_REF] Mémin | Dense estimation and object-based segmentation of the optical flow with robust techniques[END_REF]. Adapted optimization schemes must then be adopted to cope with non-linearity or non-convexity induced by the robust terms, as will be discussed in Sections 4.2 and 5.2. A priori smoothness assumption based on parametric constraint (Section 4) or explicit regularization (Section 5) also counterbalances local invalidity of data constancy.

Robust statistics and regularization treat the problem of violations of the constancy assumption by considering it as noise, with underlying distribution assumptions [START_REF] Krajsek | On the equivalence of variational and statistical differential motion estimation[END_REF][START_REF] Simoncelli | Probability distributions of optical flow[END_REF]. The considered distributions may not suitably model the possibly large localized violations implied by the above listed causes. Therefore, a large number of alternatives to brightness constancy have been proposed, aiming at more stable invariance properties. A few experimental studies have compared performances of different data costs given fixed optimization and regularization contexts [START_REF] Steinbrucker | Advanced data terms for variational optic flow estimation[END_REF][START_REF] Vogel | An evaluation of data costs for optical flow[END_REF].

Let us notice as a preamble to this section that it is difficult in practice to design a data term independently from the spatial coherence constraint and the optimization strategy to which it will be associated. For example, sophisticated feature conservation usually involves specific optimization difficulties, and thus requires careful design of the optimization solution. We will dedicate this section to a review of the main classes of data terms for optical flow estimation, emphasizing their validity domains and their limitations, independently from the estimation context in which they were elaborated.

Beyond brightness constancy

We explore several matching costs aiming at overcoming the drawback of the brightness constancy, in particular its sensitivity to noise and illumination changes.

Image transformation

A first class of data potentials exploits the same pixel-wise form as [START_REF] Aubert | Computing optical flow via variational techniques[END_REF], but operates on a tranformed version f (I) of the original image sequence:

ρ data (x, I 1 , I 2 , w) = φ( f (I 2 )(x + w(x)) -f (I 1 )(x)) (8)
Image smoothing. We can first notice that Gaussian smoothing is applied as a pre-processing step by most methods [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Zimmer | Optic flow in harmony[END_REF], in order to reduce the influence of noise. It can be viewed as a modified version of brightness constancy, setting f as a Gaussian filtering operator.

High-order constancy. Image derivatives possess illumination invariance properties that are well suited for motion estimation. The constancy of spatial image gradient, defined by f (I) = ∇I, has been introduced in [START_REF] Uras | A computational approach to motion perception[END_REF] for its ability to overcome the aperture problem when the determinant of the Hessian is nonzero. However, when applied on the directional derivative vectors, the gradient conservation only holds for translational or divergence motions. To achieve rotational invariance, the penalty should rather be applied on the magnitude of the derivatives, that is f (I) = ∇I [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF]. It was subsequently used in the context of the local approach (see Section 4) in [START_REF] Tistarelli | Multiple constraints to compute optical flow[END_REF] and integrated in global variational methods (see Section 5) in [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF]. Despite a demonstrated performance gain in the case of additive illumination changes compared to brightness constancy, gradient information is also much more sensitive to noise, and disappears in poorly textured regions. Therefore, it is always used in addition to brightness constancy. A large number of methods rely on this combination and achieve good results [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[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]. Finally [START_REF] Papenberg | Highly accurate optic flow computation with theoretically justified warping[END_REF] Texture. Another way to obtain robustness against illumination changes is to work with the cartoon and texture components of the image, as initiated in [START_REF] Wedel | An improved algorithm for TV-L1 optical flow[END_REF]. The decomposition proposed in [START_REF] Aujol | Structure-texture image decomposition -modeling, algorithms, and parameter selection[END_REF] consists in first obtaining the structure part I S by discontinuity-preserving smoothing (using the ROF model [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] in [START_REF] Wedel | An improved algorithm for TV-L1 optical flow[END_REF]), and then deriving texture part I T by subtracting I S to the original image. Additive illumination changes only affect the cartoon image while the texture image is less impacted. However, even if I T captures most image information, it is still insufficiently discriminative in some regions, and it is also more sensitive to noise than I S . To limit this drawback, the texture image used to compute optical flow is usually blended with the cartoon part by a parameter γ [START_REF] Krähenbühl | Efficient nonlocal regularization for optical flow[END_REF][START_REF] Sun | Secrets of optical flow estimation and their principles[END_REF][START_REF] Wedel | Structure-and motion-adaptive regularization for high accuracy optic flow[END_REF]:

f (I) = I T -γI S .
Colour spaces. When dealing with colour images, several photometric-invariant colour spaces can be exploited. In particular, multiplicative illumination invariance is essential for realistic illumination models [START_REF] Van De Weijer | Robust optical flow from photometric invariants[END_REF] and is achieved in the HSV space by the hue channel (local and global changes) and the saturation channel (only global changes) [START_REF] Mileva | Illumination-robust variational optical flow with photometric invariants[END_REF]. As for previously mentioned image transformations, the benefit in illumination change regions coincide with a loss of information in other parts, and the colour channels are in practice combined with the intensity valued channel [START_REF] Zimmer | Optic flow in harmony[END_REF]. Other colour spaces like normalized RGB [START_REF] Golland | Motion from color[END_REF] or spherical space [START_REF] Van De Weijer | Robust optical flow from photometric invariants[END_REF] have been investigated.

Combination of linear filters. As mentioned before, it is often necessary to combine several constancy assumptions. It has been achieved in [START_REF] Weber | Robust computation of optical flow in a multiscale differential framework[END_REF] by considering K different linear filters f k providing a system of constancy equations to be solved at each pixel. In [START_REF] Sun | Learning optical flow[END_REF], a learning approach is proposed for finding the best combination of constraints. The data constraints induced by the linear filters f k are imposed through a Gaussian Scale Mixture (GSM) φ GS M . [START_REF] Ayer | Layered representation of motion video using robust maximum-likelihood estimation of mixture models and MDL encoding[END_REF] where * denotes the convolution operator and Ξ k are the parameters of the GSM associated to each filter. By learning the parameters of the GSM from a set of ground truth sequences, the "weights" of each filter in (9) are automatically determined. These weights are constant in the whole image. Spatially adaptive combination is also of upmost importance and will be addressed in Section 3.2.1.

ρ data (x, I 1 , I 2 , w) = K k=1 φ GS M ( f k * I 2 (x + w(x)) -f k * I 1 (x), Ξ k ),

Patch-based similarity measures

Rather than pre-filtering images, neighborhood information can be integrated directly in the data term by patch-based similarity measures. Let us point out that a major issue with patchbased measure is the determination of the size an possibly shape of the patch.

Filtering the data term. In addition to pre-smoothing the images (8), Bruhn et al. [START_REF] Bruhn | Lucas/Kanade meets horn/schunck: Combining local and global optic flow methods[END_REF] proposed to filter the data potential as follows:

ρ(x, I 1 , I 2 , w) = f (φ(I 2 (x + w(x)) -I 1 (x))) , ( 10
)
where f is chosen as a Gaussian filter. While it was demonstrated beneficial for very noisy sequences, this choice also significantly blurs motion discontinuities and degrades the overall performance for low amount of noise, compared to pixel-wise data term, as emphasized in [START_REF] Zimmer | Optic flow in harmony[END_REF]. This limitation is addressed in [START_REF] Drulea | Total variation regularization of localglobal optical flow[END_REF][START_REF] Rashwan | Variational optical flow estimation based on stick tensor voting[END_REF] by replacing the Gaussian filtering with anisotropic discontinuity-preserving filtering (e.g. bilateral filtering in [START_REF] Drulea | Total variation regularization of localglobal optical flow[END_REF] and tensor voting in [START_REF] Rashwan | Variational optical flow estimation based on stick tensor voting[END_REF]). The relation between adaptive filtering and non-linear diffusion has also been exploited in [START_REF] Brox | Nonlinear structure tensors[END_REF].

Correlation-based measures. Similarity measures based on cross-correlation have been extensively used for various correspondence problems [START_REF] Brown | A survey of image registration techniques[END_REF]. Normalized Cross Correlation (NCC) is usually preferred for its invariance to linear illumination changes. The NCC for a window W(x) centered at pixel x is defined as

NCC(x, I 1 , I 2 , w) = y∈W(x) (I 2 (y + w(x)) -µ 2 (x + w(x)))(I 1 (y) -µ 1 (x)) υ 1 (x)υ 2 (x + w(x)) (11)
where for i = {1, 2}, µ i (x) is the mean and υ i (x) the standard deviation of I i in the window W(x). The associated data potential is

ρ data (x, I 1 , I 2 , w) = 1 -NCC(x, I 1 , I 2 , w). (12)
NCC is actually discriminative enough to be used in a matching procedure without additional regularization, and produces coarse but reasonably robust motion fields. It is used in several applications like stereovision [START_REF] Delon | Small baseline stereovision[END_REF], fluid flow analysis [START_REF] Becker | Variational adaptive correlation method for flow estimation[END_REF] or biological imaging [START_REF] Kolin | Advances in image correlation spectroscopy: measuring number densities, aggregation states, and dynamics of fluorescently labeled macromolecules in cells[END_REF] where it also enables direct physical measures for diffusion processes. The computational cost of ( 11) is a major limitation. Unlike simple cross-correlation which can be efficiently computed with Fast Fourier Transform (FFT), the computation of NCC for matching purpose cannot be easily performed in the frequency domain. In [START_REF] Lewis | Fast normalized cross-correlation[END_REF] only the numerator is computed with FFT and the denominator is rewritten as a product of sums, independent of the position of the pixel and thus efficiently computable with integral images [START_REF] Facciolo | Integral images for block matching[END_REF]. In [START_REF] Luo | A fast normalized cross-correlation calculation method for motion estimation[END_REF] this idea is generalized and the numerator is also computed with integral images, dramatically reducing the computation time and making it invariant to the patch size.

The challenge of integrating NCC in a variational optimization scheme is to cope with its non-differentiability. Indeed, Taylor expansion on the terms containing w(x) in [START_REF] Baker | Lucas-Kanade 20 years on: A unifying framework[END_REF] still yields a highly non-linear potential. The approach of [START_REF] Molnár | Illumination-robust variational optical flow using cross-correlation[END_REF][START_REF] Werlberger | Efficient minimization of the non-local Potts model[END_REF] has been applied to NCC but is able to handle arbitrary data terms as well. The authors directly linearize the data term and compute its spatial derivatives with finite differences. Werlberger et al. [START_REF] Werlberger | Efficient minimization of the non-local Potts model[END_REF] keeps a second-order approximation to ensure the convexity of the energy, necessary in the primal-dual scheme used. Another recent technique allowing fast computation of NCC relies on the fact that NCC is actually equivalent to the Sum of Square Differences (SSD) when the images are filtered with the cheaply computed correlation transform [START_REF] Drulea | Motion estimation using the correlation transform[END_REF].

Census. Census transform [START_REF] Zabih | Non-parametric local transforms for computing visual correspondence[END_REF] recently regained interest and was promoted in [START_REF] Stein | Efficient computation of optical flow using the census transform[END_REF] for optical flow estimation [START_REF] Hafner | Why is the census transform good for robust optic flow computation?[END_REF][START_REF] Hermann | Hierarchical scan-line dynamic programming for optical flow using semi-global matching[END_REF][START_REF] Mohamed | TV-L1 optical flow estimation with image details recovering based on modified census transform[END_REF][START_REF] Müller | Illumination-robust dense optical flow using census signatures[END_REF][START_REF] Muller | Feature-and depthsupported modified total variation optical flow for 3d motion field estimation in real scenes[END_REF][START_REF] Ranftl | Pushing the limits of stereo using variational stereo estimation[END_REF][START_REF] Vogel | An evaluation of data costs for optical flow[END_REF]. The Census signature is a bit string reflecting relative value of pixels of a patch with respect to the center pixel. By discarding the absolute intensity values, only the structure of the neighborhood is encoded in the signature, which makes it robust to illumination changes. It has shown robust behaviour in outdoor scenes and vehicle driving scenarios [START_REF] Ranftl | Pushing the limits of stereo using variational stereo estimation[END_REF][START_REF] Stein | Efficient computation of optical flow using the census transform[END_REF][START_REF] Vogel | An evaluation of data costs for optical flow[END_REF]. Integrating the Census transform in variational optical flow is not trivial since it cannot be easily linearized. Solutions to remedy this problem are convex approximation [START_REF] Vogel | An evaluation of data costs for optical flow[END_REF], reformulation as a generalization of the gradient constancy conservation [START_REF] Hafner | Why is the census transform good for robust optic flow computation?[END_REF] or linearization of the data term [START_REF] Müller | Illumination-robust dense optical flow using census signatures[END_REF][START_REF] Ranftl | Pushing the limits of stereo using variational stereo estimation[END_REF] as previously mentioned for NCC [START_REF] Werlberger | Efficient minimization of the non-local Potts model[END_REF]

Spatially adaptive data cost

The validity of each of matching costs is limited to a given range of visual situations. In a single frame, regions can coexist satisfying a given constancy assumption, and violating others. One solution could be to linearly combine them to take advantage of their complementary invariance properties. Softly selecting the best constancy constraint at each pixel is usually devoted to the robust penalty function, limiting the influence of local constraint violations. However, the data term should ideally be spatially adapted.

We distinguish two classes of methods achieving the spatial adaptivity: i) optimization of the weights of a linear combination of data potentials, and ii) estimation of the spatial distribution of the errors attached to a single data potential. The normalization of the data term used in [START_REF] Schoenemann | Near real-time motion segmentation using graph cuts[END_REF][START_REF] Simoncelli | Probability distributions of optical flow[END_REF][START_REF] Zimmer | Optic flow in harmony[END_REF] could fall in this category since a spatially varying weight is applied to the data term. It is derived from the linearized brightness constancy equation to prevent too strong constraints in regions of high image gradient (see a detailed interpretation in [START_REF] Zimmer | Optic flow in harmony[END_REF]).

Combination of constancy assumptions

The combination of P data constraints can be expressed as the weighted sum of their associated potentials ρ p (x, I 1 , I 2 , w):

ρ data (x, I 1 , I 2 , w) = P p=1 λ p (x) ρ p (x, I 1 , I 2 , w). (13)
Weights λ p (x) are spatially variant and have to be optimized to locally favor different data terms.

The idea of combining several data constraints has already been explored twenty years ago in [START_REF] Heitz | Multimodal estimation of discontinuous optical flow using Markov random fields[END_REF]. In addition to the classical brightness constancy, the authors exploited a complementary sparse edge-based constraint. Weights λ p (x) are binary confidence measures derived from hypothesis testing providing evidence on each constraint.

In [START_REF] Xu | Motion detail preserving optical flow estimation[END_REF], intensity and gradient conservation are combined and their complementarity is experimentally demonstrated. The weights are defined to operate a binary selection between the two constraints and are obtained by considering a mean field approximation of ( 13), which intuitively amounts to selecting the constraint having the lower (normalized) potential. This idea has been used subsequently in [START_REF] Mozerov | Constrained optical flow estimation as a matching problem[END_REF].

The work of [START_REF] Kim | Optical flow via locally adaptive fusion of complementary data costs[END_REF] addresses the problem in its most general form [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF], allowing the combination of an arbitrary number and type of data conservation assumptions. A confidence measure for arbitrary data term is designed as an extension of the feature discriminability [START_REF] Shi | Good features to track[END_REF] to data discriminability. The confidence measures are used as local constraints on the weights λ p (x) in ( 13), and a regularization on λ p (x) is also imposed. The weights are then optimized jointly with the motion field.

Modeling intensity variations

Another way to handle errors related to the constancy constraint is to design more elaborated models of intensity variations. The goal is to explicitly model the intensity changes due to environmental effects so as to isolate more accurately motion-induced changes.

The most common approach introduces an additional variable representing deviations from constancy assumption by a function e(x, I 1 , I 2 , ξ) parametrized by the vector ξ:

ρ data (x, I 1 , I 2 , w, ξ) = φ(I 2 (x+w(x))-I 1 (x)-e(x, I 1 , I 2 , ξ)). ( 14
)
The model proposed by [START_REF] Negahdaripour | Revised definition of optical flow: Integration of radiometric and geometric cues for dynamic scene analysis[END_REF] is composed of an offset change ξ o : Ω → R, accounting, e.g., for moving shadows or highlights, and a multiplicative change ξ m : Ω → R encoding linear illumination variations. The error function can then be expressed as

e(x, I 1 , I 2 , ξ) = ξ m (x)I 1 (x) + ξ o (x). ( 15
)
This general formulation has been exploited in a number of works, considering either the offset parameter alone [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Fortun | Aggregation of patch-based estimations for illumination-invariant optical flow in live cell imaging[END_REF][START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF], the multiplier alone [START_REF] Zach | Fast gain-adaptive KLT tracking on the GPU[END_REF] or both parameters [START_REF] Kim | Robust motion estimation under varying illumination[END_REF][START_REF] Lai | Robust image matching under partial occlusion and spatially varying illumination change[END_REF][START_REF] Teng | Accurate optical flow computation under non-uniform brightness variations[END_REF]. They may differ on the type of spatial coherency, the penalty function, or the optimization strategy. A smoothness constraint is assumed on ξ o and ξ m , either with a local parametric form [START_REF] Negahdaripour | Revised definition of optical flow: Integration of radiometric and geometric cues for dynamic scene analysis[END_REF][START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF] or a global regularization [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Fortun | Aggregation of patch-based estimations for illumination-invariant optical flow in live cell imaging[END_REF][START_REF] Kim | Robust motion estimation under varying illumination[END_REF][START_REF] Lai | Robust image matching under partial occlusion and spatially varying illumination change[END_REF][START_REF] Teng | Accurate optical flow computation under non-uniform brightness variations[END_REF]. The offset formulation is also used in [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], but is associated to a sparsity constraint aiming at retrieving violations due to occlusions. The more general approach of [START_REF] Demetz | Learning brightness transfer functions for the joint recovery of illumination changes and optical flow[END_REF] parametrizes intensity changes in terms of Brightness Transfer Function [START_REF] Grossberg | Modeling the space of camera response functions[END_REF], which coefficients are learned from training data.

The model ( 15) is based on a general polynomial approximation. If specific knowledge about the observed physical process is available, dedicated models can be designed. A number of physical constraints are explored in [START_REF] Haussecker | Computing optical flow with physical models of brightness variation[END_REF], where a generic local estimation framework is designed based on a Taylor expansion of arbitrary data constraints similar to the subsequent methods [START_REF] Molnár | Illumination-robust variational optical flow using cross-correlation[END_REF][START_REF] Werlberger | Efficient minimization of the non-local Potts model[END_REF].

Parametric approach

As explained in Section 2.4, a spatial coherence constraint on the flow field must be added to the previously described data terms. To this end, one can impose the flow field to follow a parametric model in region R ⊆ Ω. The motion field w α : R → R 2 is then fully characterized by the associated parameter vector α. When the region R is a small sub-domain of the image, these methods are referred to as local approaches. The objective energy to be minimized is the weighted sum of the potentials provided by each pixel of R:

α = arg min α x∈R g(x) ρ data (x, I 1 , I 2 , w α ) (16)
where g(x) is a spatial weighting function controlling the influence of pixel x in the estimation.

It is crucial to determine the local estimation domain where the parametric form of the motion model is a valid approximation of the true motion. Low-order polynomial motion models like translation or affine deformation can usually represent motion in small neighborhoods, whereas more complex models like deviations from affine constraint or combination of basis functions can deal with larger regions.

We first give an overview of the mostly used motion models and their associated optimization strategies. Secondly, we discuss about the different ways to define appropriate local estimation supports.

Motion models

The choice of the motion model is driven by a trade-off between efficiency and representativeness. Complex nonlinear and physical-based models can be exploited to model deformations for image registration [START_REF] Sotiras | Deformable medical image registration: A survey[END_REF]. These models are particularly well adapted to physically constrained situations as they can be encountered in medical imaging, and in particular to capture smooth deformations. In contrast, optical flow is dealing with temporal sequences of arbitrary content, usually involving motions of several objects with unrelated behaviours, generating motion discontinuities as well as smooth motion parts. As a result, it is difficult to capture the whole complexity of motion fields with a single unifying and computationally tractable parametric model. Therefore, attempts in this direction are not frequent and not among the best performing methods in optical flow benchmarks. The approach of most parametric methods for optical flow is rather to rely on much simpler motion models, mostly polynomial models, and to restrict their application to local domains, where they can represent accurate approximations.

We will restrict ourselves to linear models of the form

w α (x) = K k=1 b k (x)α k (x), ( 17
)
where b = {b k (x)} k∈[1,K] , b k (x) ∈ R 2 are basis functions and α(x) = {α k (x)} k∈[1,K]
are the parameters to be optimized. Other parametric models than those described here can be found, like planar surfaces or rigid body [START_REF] Bergen | Hierarchical model-based motion estimation[END_REF], or wavelet bases [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF][START_REF] Dérian | Wavelets and optical flow motion estimation[END_REF][START_REF] Shen | Sparsity model for robust optical flow estimation at motion discontinuities[END_REF][START_REF] Wu | Image registration using wavelet-based motion model[END_REF]. It can also be noted that parametric models are sometimes completed with explicit regularization terms (see Section 5) imposed on the parameters themselves [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF][START_REF] Ju | Skin and bones: Multi-layer, locally affine, optical flow and regularization with transparency[END_REF][START_REF] Mémin | Hierarchical estimation and segmentation of dense motion fields[END_REF][START_REF] Nir | Over-parameterized variational optical flow[END_REF].

Polynomial models

Polynomial models are among the most compact parametric representations of motion fields and are also remarkably well suited to retrieve local physical motion of individual objects and even deformable motion. The polynomial model can be seen as a special case of ( 17), which we write for clarity as:

w θ (x) = B(x) θ, (18)
where B(x) is a matrix determining the form of the model and θ is the parameter vector. Apart from the exception of [START_REF] Nir | Over-parameterized variational optical flow[END_REF] where the parameters are spatially variant and regularized, the parameter θ is kept constant over the estimation domain. The matrix B(x) depends on the pixel coordinates x ∈ R, and the order of the polynomial determines the complexity of the motion field. Low-order polynomials are usually sufficient to model smooth motion fields, and their small number of parameters allows for efficient computation. The two mostly used polynomial models correspond to the translation and the affine motion:

Translational : θ = (a 1 , a 2 ) ; B(x) = 1 0 0 1 Affine : θ = (a 1 , a 2 , a 3 , a 4 , a 5 , a 6 ) ; B(x) = 1 x 1 x 2 0 0 0 0 0 0 1 x 1 x 2 with x = (x 1 , x 2 ) .
The translation assumption is very restrictive and must be applied to very small regions [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF]. The physical assumption underlying the affine model is a rigid motion of 3D planar surfaces projected orthogonally on the image plane. Higher-order polynomials can model more complex situations, but are still too smooth to allow for motion discontinuities. For example, the 8-parameter quadratic model represents rigid motion of a planar surface in perspective projection. The small number of parameters of the affine model and its realistic local assumption make it often considered as the best trade-off between complexity and descriptiveness [START_REF] Black | The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields[END_REF][START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF][START_REF] Mémin | Hierarchical estimation and segmentation of dense motion fields[END_REF][START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF]. The accuracy of of local affine motion estimation in appropriate estimation support was experimentally demonstrated in [START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF].

Learned basis

The basis functions can also be learned from a set of training flow fields. As for polynomial models, the resulting learned motion models cannot describe motion of complex scenes, but they are able to retrieve a larger diversity of local motion patterns, including discontinuities.

The design of the training set reflects the assumption on the form of the flow. In a generic point of view, Black et al. [START_REF] Black | Learning parameterized models of image motion[END_REF] used synthetic motion fields representing simple motion patterns. The method of [START_REF] Nieuwenhuis | Complex motion models for simple optical flow estimation[END_REF] relies on a large number of patches of ground truth motion fields. The training set can be dedicated to a specific application, as in [START_REF] Fleet | Design and use of linear models for image motion analysis[END_REF] where the aim is to estimate mouth motion. To avoid resorting to external ground truth, [START_REF] Garg | Robust trajectory-space TV-L1 optical flow for non-rigid sequences[END_REF] defines the training set on the processed sequence itself. The basis set is composed of trajectories constructed by feature tracking on large temporal scales, in regions ensuring reliable tracks. In all these works, an orthogonal basis of flow fields is generated by PCA decomposition, conserving only the first K components containing most of the variance of the training set.

Free-form deformations

The free-form deformation model (FFD) [START_REF] Rueckert | Nonrigid registration using free-form deformations: application to breast MR images[END_REF] has been originally introduced for image registration and has demonstrated great efficiency to retrieve smooth deformations. The displacements are defined on a coarse regular subgrid of the image and is interpolated on the final resolution with B-splines. The motion basis b k (x) is thus formed with the displacements of the K control points and coefficients α k (x) are B-spline influence functions. The dimensionality reduction induced by the subsampling of the image grid makes the computation much easier, and the spatial coherence of the deformation is ensured by the B-spline interpolation. On the counterpart, the framework cannot retrieve sharp motion discontinuities, while it is necessary for optical flow applications. Image-adaptive non-regular control points distribution [START_REF] Schnabel | A generic framework for non-rigid registration based on nonuniform multi-level free-form deformations[END_REF] or coarse-to-fine spacing strategies [START_REF] Rueckert | Nonrigid registration using free-form deformations: application to breast MR images[END_REF] are possibilities to address this issue.

The work of [START_REF] Szeliski | Motion estimation with quadtree splines[END_REF] was the first to apply this idea to optical flow, with non-uniform control points defined on imagedriven quadtrees. In [START_REF] Pingault | A robust multisacle B-spline function decomposition for estimating motion transparency[END_REF], B-Splines defined on a uniform grid are used to retrieve smooth deformations. In [START_REF] Glocker | Optical flow estimation with uncertainties through dynamic MRFs[END_REF] the problem is turned in a discrete setting. The range of motion labels is iteratively adapted by estimating a local uncertainty covariance. Despite good results, the method is still limited by oversmoothing. The method of [START_REF] Shi | Registration using sparse free-form deformations[END_REF] addresses the discontinuity problem with a sparsity constraint on the B-spline coefficients, allowing to modulate the influence of the control points. Recent approaches [START_REF] Glocker | Inter and intra-modal deformable registration: Continuous deformations meet efficient optimal linear programming[END_REF][START_REF] Glocker | Optical flow estimation with uncertainties through dynamic MRFs[END_REF][START_REF] Shi | Registration using sparse free-form deformations[END_REF] also add an explicit regularization on the motion field, overlapping with the methodology described in Section 5.

Optimization

Parametric models are usually associated with the penalty of a pixel-wise data constancy constraint [START_REF] Aubert | Computing optical flow via variational techniques[END_REF]. In case of intensity constancy [START_REF] Amiaz | Coarse to over-fine optical flow estimation[END_REF], energy [START_REF] Barnes | The generalized patchmatch correspondence algorithm[END_REF] then writes:

E(α) = x∈R g(x)φ (∇I(x) • w α (x) + I t (x)) ( 19
)
where ∇I is the image gradient and I t = ∂I/∂t. The special case of a quadratic function φ and a translational model as in [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF] leads to a very simple optimization problem, since the cancelling of the derivatives of ( 19) amounts to solving the linear system Mα = b, with

M = x g(x) I 2 x 1 (x) x g(x) I x 1 (x)I x 2 (x) x g(x) I x 1 (x)I x 2 (x) x g(x) I 2 x 2 (x) b = - x g(x) I x 1 (x)I t (x) x g(x) I x 2 (x)I t (x)
where I x 1 and I x 2 are the partial derivatives of I respectively along the horizontal axis x 1 and the vertical axis x 2 . The rank of M allows one to decide if a unique solution of the linear system exists, and can be used to adapt the size of the local domain R (see Section 4.3.2). Despite the limitations of the quadratic penalty, this approach has become very popular for its implementation simplicity, low computational cost and available code in the OpenCV library [START_REF] Bouguet | Pyramidal implementation of the Lucas-Kanade feature tracker. Description of the algorithm[END_REF]41]. However, robust estimation is often advocated [START_REF] Black | The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields[END_REF][START_REF] Black | Learning parameterized models of image motion[END_REF][START_REF] Dérian | Wavelets and optical flow motion estimation[END_REF][START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF][START_REF] Senst | Robust local optical flow for feature tracking[END_REF] as mentioned in Section 3, especially for polynomial models, to deal with the frequent case of multiple motions in the estimation domain. Among the variety of optimization methods used in case of robust penalty function, the Iterative Reweighted Least Squares (IRLS) [START_REF] Holland | Robust regression using iteratively reweighted least-squares[END_REF] and gradient descent approaches have mostly been used. IRLS proceeds by successive optimizations of quadratic problems weighted by a function of the current estimate and is implemented in the Motion2D software [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF][START_REF] Senst | Robust local optical flow for feature tracking[END_REF] with . Gradient descent approaches are often coupled with Graduated Non-Convexity (GNC) [START_REF] Black | The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields[END_REF][START_REF] Blake | Visual reconstruction[END_REF] to cope with non-convexity of [START_REF] Basset | Recovery of motion patterns and dominant paths in videos of crowded scenes[END_REF]. Regarding the slow convergence of steepest descent, it is preferable to use second order approximations and Newton methods, or quasi Newton methods like L-BFGS or Levenberg-Marquardt, approximating the Hessian for large dimension problems.

Neighborhood selection

As previously mentioned, the spatial adaptation of the parameters α(x) is a way to cope with complex and discontinuous flow fields [START_REF] Garg | Robust trajectory-space TV-L1 optical flow for non-rigid sequences[END_REF][START_REF] Nir | Over-parameterized variational optical flow[END_REF][START_REF] Shi | Registration using sparse free-form deformations[END_REF]. This approach often involves a priori constraints on the spatial distribution of α(x) and is thus strongly related to the methods that will be presented in Section 5. Such a dense parameter map moves away from the compactness of parametric models.

On the other hand, when parameters are constant over the estimation domain R, the resulting motion field is smooth and can constitute a valid approximation only in regions of coherent motion. The region R must be large enough to enable motion estimation, while small enough to keep valid the parametric approximation (generalized aperture problem [START_REF] Black | The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields[END_REF]). We will describe below the strategies for defining R in the case of constant parameters α(x) = α over R, for polynomial models.

Entire domain

Despite their inability to retrieve realistic motion in arbitrary scenes, polynomial models combined with robust estimation are well adapted to capture the dominant image motion. Applied in the whole image domain Ω, they become particularly useful to separate or compensate the camera motion [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF]. Applications to moving object detection [START_REF] Odobez | Separation of moving regions from background in an image sequence acquired with a mobile camera[END_REF] or to action recognition [START_REF] Piriou | Recognition of dynamic video contents with global probabilistic models of visual motion[END_REF] are among the most common ones.

Square patches

The approach initiated by [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF] performs independent estimations in small square or circular patches. Most of the related methods use fixed patch size, and conserve the velocity vector deduced from the estimated motion model at the square center [START_REF] Baker | Lucas-Kanade 20 years on: A unifying framework[END_REF][START_REF] Bigun | Multidimensional orientation estimation with applications to texture analysis and optical flow[END_REF][START_REF] Kim | A local approach for robust optical flow estimation under varying illumination[END_REF][START_REF] Zach | Fast gain-adaptive KLT tracking on the GPU[END_REF]. This choice is very popular for its simplicity of implementation and can be naturally parallelized [START_REF] Sinha | Feature tracking and matching in video using programmable graphics hardware[END_REF][START_REF] Zach | Fast gain-adaptive KLT tracking on the GPU[END_REF]. It is still extensively used for numerous applications. However, following the generalized aperture problem, patches centered at each pixel with a fixed size are likely to contain either multiple motions or no image gradient.

Multiple motions in a single patch can be partially handled with robust estimation by rejecting secondary motions, considered as outliers [START_REF] Black | The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields[END_REF][START_REF] Gelgon | A region tracking method with failure detection for an interactive video indexing environment[END_REF][START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF][START_REF] Senst | Robust local optical flow for feature tracking[END_REF].

The second option is to adapt the size or the position of the patch so that it contains an unimodal motion distribution. The size of the patch can be adapted with a bias-variance criterion [START_REF] Maurizot | Determination of singular points in 2d deformable flow fields[END_REF], or based on a confidence measure on the reliability of the local domain for parametric estimation. Starting from a small patch size, it is thus possible to increase the size of the patch until the condition for a reliable domain is violated [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]. In the translational model case, one can analyze the singularity of the matrix M of the linear system (20) and, for example, impose a minimum threshold to the maximum eigenvalue of M [START_REF] Barron | Evaluation of optical flow[END_REF]. Rather than adapting the size, [START_REF] Jodoin | Optical-flow based on an edge-avoidance procedure[END_REF] adapts the position of the patches. Patches corrupted by strong intensity edges are displaced by a mean-shift procedure to reach more homogeneous regions. The method described in [START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF] adapts both size and position of the square patches. This adaptation is performed thanks to the aggregation framework decoupling the estimation of candidates with predefined sizes and positions, and their selection achieved with a global model (see Section 5).

Finally, local adaptation of the weights g(x) in ( 16) is another way to adapt the estimation support. It can be done for instance with bilateral filter weights as in [START_REF] Drulea | Total variation regularization of localglobal optical flow[END_REF], or using the non-linear structure tensor [START_REF] Brox | Nonlinear structure tensors[END_REF].

Segmented regions

The optimal supports to perform polynomial motion estimation ideally correspond to a segmentation of the image in coherently moving regions. We briefly describe two types of approaches: independent image segmentation and joint estimation of motion and region supports or frontiers.

Image segmentation. While the ultimate goal is to segment the unknown motion field, color-based image segmentation is a much simpler alternative which can help motion estimation. It can be reasonably assumed that motion discontinuities coincide with image discontinuities (but the inverse is far from being true). It implies that an image segmentation is a motion field over-segmentation, and obtained regions are thus guaranteed to contain no motion discontinuity. However, merely estimating motion in the resulting regions is problematic for two reasons.

The first limitation is that the segmented regions may not contain enough information for motion estimation. Parametric estimations in these regions must be performed by circumvented ways. The very fine over-segmentation of [START_REF] Zitnick | Consistent segmentation for optical flow estimation[END_REF] imposes for instance to perform region matching. Generally, an independent coarse and cheap motion estimation is fused with the color image segmentation to overcome the lack of information [START_REF] Black | Estimating optical flow in segmented images using variable-order parametric models with local deformations[END_REF][START_REF] Bleyer | Segmentation-based motion with occlusions using graph-cut optimization[END_REF][START_REF] Xu | A segmentation-based variational model for accurate optical flow estimation[END_REF]. In [START_REF] Xu | A segmentation-based variational model for accurate optical flow estimation[END_REF], hybrid regions are found by applying mean-shift segmentation in the extended space of color and motion. Differently, [START_REF] Black | Estimating optical flow in segmented images using variable-order parametric models with local deformations[END_REF][START_REF] Bleyer | Segmentation-based motion with occlusions using graph-cut optimization[END_REF] fit a parametric flow field on the coarse initial motion field, obtained with a global regularized method [START_REF] Black | A framework for the robust estimation of optical flow[END_REF] for [START_REF] Black | Estimating optical flow in segmented images using variable-order parametric models with local deformations[END_REF], and with the sparse KLT tracker [START_REF] Shi | Good features to track[END_REF] for [START_REF] Bleyer | Segmentation-based motion with occlusions using graph-cut optimization[END_REF].

The second problem is that spatial coherence between estimated motion in neighboring segments is not ensured. Global regularization (see Section 5) can here be imposed, either on the motion parameters associated to each region [START_REF] Xu | A segmentation-based variational model for accurate optical flow estimation[END_REF] (similarly to [START_REF] Ju | Skin and bones: Multi-layer, locally affine, optical flow and regularization with transparency[END_REF], not resorting to image segmentation), on the coarse motion field completing color information [START_REF] Black | Estimating optical flow in segmented images using variable-order parametric models with local deformations[END_REF] or, in a layered approach, on the layer assignment function [START_REF] Bleyer | Segmentation-based motion with occlusions using graph-cut optimization[END_REF] Joint estimation and segmentation. Color segmentation is usually too dependent on the image content to make it the basis of a robust motion estimation method. Rather than considering segmentation and estimation as two independent tasks, most methods have a coupled approach of the problem. Motion parameters and region supports are jointly estimated by minimizing a global energy imposing a coupling between them. This approach has first been addressed as a labelling problem [START_REF] Bouthemy | Motion segmentation and qualitative dynamic scene analysis from an image sequence[END_REF][START_REF] Odobez | Direct incremental model-based image motion segmentation for video analysis[END_REF] where the label field l : Ω → {l 1 , . . . , l N } associated to the N regions is estimated jointly with the motion parameters in each region α = {α 1 , . . . , α N }, in a discrete Markov Random Field framework:

E(α, l) = x∈Ω d ρ data x, I 1 , I 2 , w α l(x) + <x,y> ρ MRF reg (l(x), l(y)), (20)
where Ω d is the discrete image domain and ρ MRF reg (l(x), l(y)) is a regularization prior on the label field, typically chosen as ρ MRF reg (l(x), l(y)) = 1 -δ(l(x), l(y)), with δ(•) the Kronecker function. The optimization is done alternatively between motion and regions. Another viewpoint in a variational framework extends the Mumford-Shah formulation of image segmentation [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] to motion segmentation [START_REF] Cremers | Motion competition: A variational approach to piecewise parametric motion segmentation[END_REF][START_REF] Paragios | Geodesic active regions and level set methods for motion estimation and tracking[END_REF][START_REF] Vásquez | Joint multiresolution segmentation and parametric estimation of image motion by basis function representation and level set evolution[END_REF]. In addition to the data fitting potential inside each region R i , a constraint restricting the length L(C) of the set of region boundaries C is imposed globally, resulting in the energy

E(α, C) = N i=1 R i ρ data (x, I 1 , I 2 , w α i )dx + ν L(C). ( 21
)
The minimization is performed alternatively on the flow and the boundaries. Minimizing [START_REF] Beauchemin | The computation of optical flow[END_REF] with respect to C requires a differentiable approximation of the contour length L(C). It is common to implicitly represent the partitioning of the image with level sets, which allows to represent the interior of the regions by the sign of the function, as well as the total length of the boundaries by their level lines. One level set function can only represent two regions. For an arbitrary number of N regions, it is possible to define N corresponding levels sets, at the price of a high computational cost and a more complex energy to prevent vacuums in the partitioning. Other strategies can be employed, as the one of [START_REF] Chan | Euler's elastica and curvaturebased inpainting[END_REF] re-used in [START_REF] Cremers | Motion competition: A variational approach to piecewise parametric motion segmentation[END_REF], for more sophisticated combinations between functions. In [START_REF] Paragios | Geodesic active regions and level set methods for motion estimation and tracking[END_REF], this level set framework is augmented with an edge-driven tracking and background detection. In [START_REF] Riklin-Raviv | Shape-based mutual segmentation[END_REF], the two images are jointly segmented, and they influence each other through a dynamic prior term determined by the similarity of the two evolving shapes. Segmentation of static and dynamic textures is performed in [START_REF] Fazekas | Dynamic texture detection based on motion analysis[END_REF] by assigning different data conservation assumptions to each class. A graph-cut optimization scheme has been proposed in [START_REF] Dupont | Extraction of layers of similar motion through combinatorial techniques[END_REF]. These two formulations ( 20) and ( 21) can be found in numerous other works [START_REF] Kervrann | Multiscale neighborhood-wise decision fusion for redundancy detection in image pairs[END_REF][START_REF] Mémin | Hierarchical estimation and segmentation of dense motion fields[END_REF][START_REF] Schnörr | Motion-based identification of deformable templates[END_REF][START_REF] Unger | Joint motion estimation and segmentation of complex scenes with label costs and occlusion modeling[END_REF]. Similarly, layered approaches [START_REF] Ayer | Layered representation of motion video using robust maximum-likelihood estimation of mixture models and MDL encoding[END_REF][START_REF] Smith | Layered motion segmentation and depth ordering by tracking edges[END_REF][START_REF] Sun | Secrets of optical flow estimation and their principles[END_REF][START_REF] Sun | Layered segmentation and optical flow estimation over time[END_REF][START_REF] Wang | Representing moving images with layers[END_REF][START_REF] Xiao | Motion layer extraction in the presence of occlusions using graph cuts[END_REF] aim at decomposing the observed scene into overlapping layers and introduce information about depth ordering of regions. An important interest of such methods is to provide a natural occlusion handling framework.

Two main drawbacks affect the joint motion estimation and segmentation approach. First, alternating minimization w.r.t. regions and flow fields is computationally expensive. The related layered approach [START_REF] Sun | Layered image motion with explicit occlusions, temporal consistency, and depth ordering[END_REF][START_REF] Sun | Layered segmentation and optical flow estimation over time[END_REF] achieving state-of-the-art results requires several hours to process a pair of 640×480 pixels, and even GPU-based implementation [START_REF] Unger | Joint motion estimation and segmentation of complex scenes with label costs and occlusion modeling[END_REF] can need up to an hour.

Second, the initialization of motion and segmentation parameters usually has a substantial impact on the evolution of the region contours throughout the minimization procedure. Therefore, paradoxically, the best performing methods in terms of accuracy of optical flow fields resort to pre-calculated motion fields obtained with independent methods to initialize their algorithm. For example, the results reported in [START_REF] Sun | Layered image motion with explicit occlusions, temporal consistency, and depth ordering[END_REF] and [START_REF] Unger | Joint motion estimation and segmentation of complex scenes with label costs and occlusion modeling[END_REF] are obtained with an initialization by the result of the methods of [START_REF] Sun | Secrets of optical flow estimation and their principles[END_REF] and [START_REF] Werlberger | Anisotropic Huber-L1 optical flow[END_REF] respectively. Nevertheless it is shown in [START_REF] Unger | Joint motion estimation and segmentation of complex scenes with label costs and occlusion modeling[END_REF] that the initialization affects more the speed of convergence than the final result. Some works exploit more complex flow field representations than polynomial approximation, by authorizing deformations from an initial affine model [START_REF] Black | Estimating optical flow in segmented images using variable-order parametric models with local deformations[END_REF][START_REF] Mémin | Hierarchical estimation and segmentation of dense motion fields[END_REF][START_REF] Sun | Layered image motion with explicit occlusions, temporal consistency, and depth ordering[END_REF] or an explicit regularized model [START_REF] Amiaz | Piecewise-smooth dense optical flow via level sets[END_REF][START_REF] Brox | Variational motion segmentation with level sets[END_REF]. Allowing such complex and discontinuous motion fields in segmented regions actually tends too produce larger regions, and ultimately leads to global approaches where motion discontinuities are handled by the global modeling itself. Consequently, regions may not represent coherent motion, but rather a delineation adapted to the specific estimation method. Finally, imposing prior on shapes and contours has been investigated in medical imaging [START_REF] Kervrann | Statistical deformable model-based segmentation of image motion[END_REF], but a limited number of regions can be handled at the same time.

Regularized models

In an alternating manner to the parametric representation of spatial coherence, mostly adapted to smooth deformations, the global form of the motion field can be imposed by an explicit regularization term. Motion discontinuities are then no more represented by the boundaries of the regions delimiting parametric motion fields, but they are involved in the global model, often considered as outliers w.r.t. smoothness assumptions. The variational approach has been initially proposed by [START_REF] Horn | Determining optical flow[END_REF] and is usually referred to as the global approach, since the regularization term interconnects all the pixels of the image and thus requires the optimization of the objective energy to be performed globally. In this subsection, we review current versions of the regularization model and optimization strategies.

In its most general form, the energy minimized by global regularization methods can be written as:

E global (w) = Ω ρ data (x, I 1 , I 2 , w) + λ ρ reg (x, w)dx (22)
where ρ data (x, I 1 , I 2 , w) is the data potential, as discussed in Section 3, ρ reg (x, w) is the regularization potential encoding an a priori assumption on the field w, and λ is a parameter tuning the balance between the two terms. Broadly speaking, the regularization potential aims at smoothing the motion field in regions of coherent motion while preserving motion discontinuities at the boundaries of moving objects. Finding the trade-off can also be partially addressed in the adaptation of the balance parameter λ [START_REF] Héas | Bayesian inference of models and hyperparameters for robust optical-flow estimation[END_REF][START_REF] Krajsek | Bayesian model selection for optical flow estimation[END_REF][START_REF] Ng | Errors-in-variables modelling in optical flow problems[END_REF][START_REF] Zimmer | Optic flow in harmony[END_REF].

A major interest of the global variational framework is its versatility, allowing one to model different forms of flow fields by combining different data and regularization terms. One must nevertheless keep in mind that minimizing [START_REF] Becker | Variational adaptive correlation method for flow estimation[END_REF] is often a tricky task. The potentially unlimited combinations of data terms and regularization terms is restricted in practice to those compatible with efficient minimization. Besides, advances in optical flow have often been correlated with new possibilities offered by optimization techniques. For example, efficient Primal-Dual minimization for Total Variation regularization [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] have motivated a number of optical flow models [START_REF] Unger | Joint motion estimation and segmentation of complex scenes with label costs and occlusion modeling[END_REF][START_REF] Werlberger | Motion estimation with nonlocal total variation regularization[END_REF][START_REF] Zach | A duality based approach for realtime TV-L1 optical flow[END_REF]. The development of efficient discrete optimization techniques based on graph cuts [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF] or message passing [START_REF] Kolmogorov | Convergent tree-reweighted message passing for energy minimization[END_REF] also inspired various works [START_REF] Chen | Large displacement optical flow from nearest neighbor fields[END_REF][START_REF] Lempitsky | Fusion moves for Markov random field optimization[END_REF][START_REF] Mozerov | Constrained optical flow estimation as a matching problem[END_REF]. Another consequence of the close intricacy between energy model and optimization technique is the difficulty to compare performances of different models, as global optimum is in general not guaranteed for sophisticated energies and the quality of the local optimum depends on the type of optimization method.

We will detail in Section 5.1 existing regularization models independently from optimization techniques, for the sake of clarity. Section 5.2 focuses on the dependency between specific energy models and optimization methods.

Regularization models 5.1.1. Spatial flow gradient constraint

The most natural and widely used way to impose smoothness of the motion field is to penalize the magnitude of the flow gradient:

ρ reg (x, w) = h(x, I 1 ) φ( ∇w(x) 2 ) ( 23
)
where φ(•) is the penalty function and h(x, I 1 ) is a weighting function.

A taxonomy of optical flow regularizers has been proposed in [START_REF] Weickert | A theoretical framework for convex regularizers in PDE-based computation of image motion[END_REF]. The authors focus on convex and rotational invariance properties, and prove uniqueness of the solution in each case. For each regularization of type [START_REF] Bergen | Hierarchical model-based motion estimation[END_REF], they show the equivalence between the resolution of the Euler-Lagrange equations associated with energy [START_REF] Becker | Variational adaptive correlation method for flow estimation[END_REF] and diffusion filtering. In addition, a diffusion tensor is derived for each particular variation of [START_REF] Bergen | Hierarchical model-based motion estimation[END_REF]. We will give a more succinct overview, taking only some elements from this classification and integrating more recent approaches.

Flow-driven regularization. In flow-driven approaches, no relation between the form of the flow field and the structure of the image is assumed. The weighting function is thus h(x, I 1 ) = 1, ∀x ∈ Ω. The seminal formulation of [START_REF] Horn | Determining optical flow[END_REF] adopts a quadratic penalization function:

ρ reg (x, w) = ∇u(x) 2 + ∇v(x) 2 , (24)
with w(x) = (u(x), v(x)) . The quadratic penalization being unable to capture motion discontinuities, the introduction of a complementary line process to handle motion discontinuities in the MRF framework [START_REF] Heitz | Multimodal estimation of discontinuous optical flow using Markov random fields[END_REF][START_REF] Konrad | Bayesian estimation of motion vector fields[END_REF], and the use of robust subquadratic penalties [START_REF] Black | The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields[END_REF][START_REF] Deriche | Optical-flow estimation while preserving its discontinuities: a variational approach[END_REF][START_REF] Mémin | Dense estimation and object-based segmentation of the optical flow with robust techniques[END_REF] have soon been employed to overcome the problem. Among the wide panel of robust functions, the popular parameter-free Total Variation (TV) prior, has interesting and useful properties [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Cohen | Nonlinear variational method for optical flow computation[END_REF][START_REF] Xu | Motion detail preserving optical flow estimation[END_REF][START_REF] Zach | A duality based approach for realtime TV-L1 optical flow[END_REF]. Contrary to most other robust norms, the TV yields a convex constraint facilitating optimization. The non-differentiability in 0 is generally alleviated by using the regularized version φ(z) = √ z 2 + 2 , where is a small constant. Associated to proximal splitting minimization, TV involves solving several ROF (Rudin-Osher-Fatemi) models [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], for which very efficient algorithms exist [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. A series of optical flow estimation methods have exploited this idea for fast and accurate minimization [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Wedel | An improved algorithm for TV-L1 optical flow[END_REF][START_REF] Werlberger | Efficient minimization of the non-local Potts model[END_REF][START_REF] Zach | A duality based approach for realtime TV-L1 optical flow[END_REF].

TV regularization actually favors piecewise constant flow fields. This framework is known to transform a smoothly varying variable to a succession of small discontinuous constant steps (staircasing artifacts). This undesirable effect can be reduced by replacing the L 1 penalization by a quadratic one for small gradient magnitude, which is the behaviour of the Huber norm [START_REF] Shulman | Regularization of discontinuous flow fields[END_REF][START_REF] Werlberger | Anisotropic Huber-L1 optical flow[END_REF]. Another possibility is to penalize higher-order derivatives of the flow, as done in [START_REF] Trobin | An unbiased secondorder prior for high-accuracy motion estimation[END_REF] for the second derivative, to favor piecewise affine flow fields. The Total Generalized Variation (TGV) [START_REF] Bredies | Total generalized variation[END_REF], generalizes L 1 penalization to arbitrarily high order derivatives. The performance gain of the second order TGV has been experimentally shown for smooth deformation conditions, for which the staircasing effect is prominent with TV regularization [START_REF] Braux-Zin | A general dense image matching framework combining direct and feature-based costs[END_REF][START_REF] Ranftl | Pushing the limits of stereo using variational stereo estimation[END_REF][START_REF] Vogel | An evaluation of data costs for optical flow[END_REF].

Despite the demonstrated performance of TV due to its algorithmic attractiveness, the real distribution of optical flow derivatives has been shown to follow a more heavy-tailed and concave distribution [START_REF] Roth | On the spatial statistics of optical flow[END_REF]. Finding good approximate solutions for non-convex priors has motivated a number of works and will be discussed in Section 5.2. When appropriate minimization strategy is available, like graph cuts or the recent work of [START_REF] Ochs | iPiano: Inertial proximal algorithm for non-convex optimization[END_REF], this kind of penalty functions has proven to yield improvements compared to the TV model. Efficient methods to minimize Potts model have also recently been investigated [START_REF] Cai | Disparity and optical flow partitioning using extended Potts priors[END_REF].

Inclusion of image gradient information. It is natural to assume a link between the motion field and its source image I 1 . As already stated in Section 4.3.3 about the relationship between motion and image segmentation, it is reasonable to consider that motion discontinuities coincide with image discontinuities delineating moving objects. This information can be incorporated in the regularization through the weighting function h(x, I 1 ) taken as a smooth decreasing function of ∇I 2 [START_REF] Alvarez | A PDE model for computing the optical flow[END_REF][START_REF] Ayvaci | Sparse occlusion detection with optical flow[END_REF][START_REF] Mozerov | Constrained optical flow estimation as a matching problem[END_REF][START_REF] Wedel | An improved algorithm for TV-L1 optical flow[END_REF][START_REF] Xu | Motion detail preserving optical flow estimation[END_REF], often defined as

h(x, I 1 ) = e -∇I(x) 2 /ς 2 ( 25
)
where ς is a parameter setting the influence of the image gradient on the regularization. Despite the risk of over-segmentation, this simple weighting strategy usually improves experimental results.

The weighting function ( 25) is isotropic since the smoothing is modulated by the same value in all the directions. This is suboptimal since we would ideally like to prevent the smoothing only across the boundaries, and allow it along them. This can be achieved by defining the regularization axes differently from the horizontal and vertical axes. The eigenvectors s 1 and s 2 of the structure tensor R ρ (x) = K ρ * [∇I(x)∇I(x) ] are well adapted since s 1 is oriented across local image edges and s 2 is orthogonal to s 1 . This idea has been first introduced in [START_REF] Nagel | An investigation of smoothness constraints for the estimation of displacement vector fields from image sequences[END_REF], which regularizer has been rewritten in [START_REF] Zimmer | Optic flow in harmony[END_REF] as:

ρ reg (x, w) = 1 ∇I(x) 2 + 2κ 2 (26) κ 2 u 2 s 1 (x) + v 2 s 1 (x) + ( ∇I(x) 2 + κ 2 ) u 2 s 2 (x) + v 2 s 2 (x)
where u s i , {i = 1, 2}, are the derivatives of u along the s i axis and κ is a regularization parameter. When κ is small, the regularization is reduced in the direction of the image gradient s 1 and strengthened along image edges s 2 depending on the image gradient magnitude ∇I(x) 2 . In [START_REF] Nagel | An investigation of smoothness constraints for the estimation of displacement vector fields from image sequences[END_REF], the eigenvectors s 1 and s 2 are obtained with a radius ρ = 0 for the Gaussian filtering of the structure tensor R ρ (x). The classical artifact produced by purely image-driven regularization is an over-fitting of the flow field on image boundaries, creating artificial motion discontinuities. To reduce the impact of image gradient in the regularization, Sun et al. proposed in [START_REF] Sun | Learning optical flow[END_REF] to keep the s 1 and s 2 directions, while suppressing the weighting on ∇I(x) 2 in (26) and employing a robust penalty function φ:

ρ reg (x, w) = φ(u 2 s 1 (x)) + φ(v 2 s 1 (x)) + φ(u 2 s 2 (x)) + φ(v 2 s 2 (x)). ( 27
)
The work of [START_REF] Zimmer | Optic flow in harmony[END_REF] proposed a generalized computation of the regularization axes, oriented to follow the data constraint rather than the image edges. In analogy with the previous approach defining s 1 , s 2 from the structure tensor, they compute the eigenvectors of a so-called regularization tensor, designed to be complementary with the data term. The approach of [START_REF] Zimmer | Optic flow in harmony[END_REF] can be generalized to data potentials built from the combination of L linear constancy constraints, that is,

ρ data (x, I 1 , I 2 , w) = L =1 φ(A (x) w(x) + B (x)). ( 28
)
For this kind of data potential, the regularisation tensor is defined by

R ρ (x) = L =1 K ρ * A (x)A (x) , (29)
where K ρ is a Gaussian kernel and s 1 ,s 2 are the eigenvectors of R ρ (x). Taking L = 1, A 1 (x) = ∇I(x) and B 1 (x) = I t (x) yields the brightness constancy constraint, and the regularization tensor reduces to the structure tensor. For more elaborated data terms, as the combination of normalized brightness and gradient constancy used in [START_REF] Zimmer | Optic flow in harmony[END_REF], the resulting axes are more consistent with the data constraints.

Non-local regularization

The gradient of the flow can only provide a local constraint on the interaction between pixels. Assuming long range interactions can capture more precisely the form of the motion field. Such non-local regularization has been recently investigated in [START_REF] Drulea | Motion estimation using the correlation transform[END_REF][START_REF] Krähenbühl | Efficient nonlocal regularization for optical flow[END_REF][START_REF] Sun | Secrets of optical flow estimation and their principles[END_REF][START_REF] Werlberger | Motion estimation with nonlocal total variation regularization[END_REF] by describing the structure of the flow in an extended neighborhood N(x) in a discrete setting as:

ρ reg (x, w) = y∈N(x) k(x, y, I 1 ) φ w(x) -w(y) 2 . ( 30
)
The weights k(x, y, I 1 ) indicate which pixel y ∈ N(x) should share a similar motion with pixel x. They are derived from the bilateral filter, favoring small distances in the spatial and color spaces [START_REF] Yoon | Adaptive support-weight approach for correspondence search[END_REF]:

k(x, y, I 1 ) = exp - x -y 2 σ 2 s - I 1 (x) -I 1 (y) 2 σ 2 c ( 31
)
where σ s and σ c control the influence of spatial distance and color similarity. This approach is image-driven in a similar way to local weighting [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF], in the sense that the smoothness is weighted by the image edges. Nevertheless, the integration on a large neighborhood mitigates the influence of local gradients and more globally exhibits the structure of the objects. It is implemented as an alternate weighted median filtering in [START_REF] Sun | Secrets of optical flow estimation and their principles[END_REF] and interpreted as a low-level soft segmentation in [START_REF] Werlberger | Efficient minimization of the non-local Potts model[END_REF].

The high-order regularization causes severe optimization difficulties discussed in Section 5.2, in particular in terms of computational cost, increasing with the size of the neighborhood N(x).

Temporal coherence

A natural idea is to extend the spatial regularization described above to the temporal dimension, assuming that motion varies smoothly across consecutive frames. Similarly with the spatial case, smoothness on the time axis can be achieved either locally, based on the temporal gradient, or taking into account a longer interval by working on trajectories.

Constraint on the temporal flow gradient. The most straightforward way to model temporal smoothness is to penalize the temporal flow gradient, analogously with the spatial flow gradient in Section 5.1.1 [START_REF] Murray | Scene segmentation from visual motion using global optimization[END_REF]. This simple extension is however unrealistic since motion of objects necessarily implies a temporal change of flow fields. Thus, regularization on the spatio-temporal direction of the motion field is more appropriate, and is achieved in [START_REF] Chin | Probabilistic and sequential computation of optical flow using temporal coherence[END_REF][START_REF] Nagel | Extending the oriented smoothness constraint into the temporal domain and the estimation of derivatives of optical flow[END_REF][START_REF] Weickert | Variational optic flow computation with a spatio-temporal smoothness constraint[END_REF][START_REF] Zimmer | Optic flow in harmony[END_REF] by extension of the spatial gradient penalties described in Section 5.1.1 to the spatio-temporal dimension. However, the performance of local temporal regularization is deceiving in most cases.

Constraint on the trajectory. Regularization based on temporal derivatives is unable to model large displacements. In this case, temporal coherence is more adequately modeled along the trajectories of objects. It was done in [START_REF] Black | Recursive non-linear estimation of discontinuous flow fields[END_REF], who does not model explicitly the trajectories but estimates temporal changes on warped images to alleviate the problem of large displacements. However, the warping is done sequentially in the forward direction and is thus prone to propagate errors. In the same spirit, Volz et al. [START_REF] Volz | Modeling temporal coherence for optical flow[END_REF] considers the same coordinate system for groups of five frames, which implies an implicit natural registration. The estimation is done jointly in all frames of the sequence, which overcomes the lack of feedback with previous frames of [START_REF] Black | Recursive non-linear estimation of discontinuous flow fields[END_REF]. The trajectory constraint in [START_REF] Garg | Robust trajectory-space TV-L1 optical flow for non-rigid sequences[END_REF] is explicitly imposed by modeling the flow field as linear combination of long-term trajectory bases, obtained from few reliable tracks.

Optimization

As previously mentioned the optimization strategy employed to minimize [START_REF] Becker | Variational adaptive correlation method for flow estimation[END_REF] has a decisive influence on the final result. We give an overview of the main continuous and discrete optimization methods and point out their adaptability to specific energy terms.

Continuous methods

Resolution of the Euler-Lagrange equations. The Euler-Lagrange equations give necessary conditions for minimizing energy of the form

E(w) = Ω F(x, w, ∇w) dx, (32)
which is the case of ( 22) with F(x, w, ∇w) = ρ data (x, I 1 , I 2 , w) + λ ρ reg (x, ∇w) when the regularization is a function of the flow gradient [START_REF] Bergen | Hierarchical model-based motion estimation[END_REF]. They provide the following system of partial differential equations:

       ∂ρ data (x,I 1 ,I 2 ,w) ∂u -λ div ∂ρ reg (x,∇w) ∂∇u = 0 ∂ρ data (x,I 1 ,I 2 ,w) ∂v -λ div ∂ρ reg (x,∇w) ∂∇v = 0, ( 33
)
where ∂• ∂z is the partial derivative operator with respect to z, div(•) is the divergence operator. The system (33) can be rewritten by introducing the diffusion tensor D accounting for the relation with diffusion equations [START_REF] Weickert | A theoretical framework for convex regularizers in PDE-based computation of image motion[END_REF]:

       ∂ρ data (x,I 1 ,I 2 ,w) ∂u -λ div (D(x) ∇u(x)) = 0 ∂ρ data (x,I 1 ,I 2 ,w) ∂v -λ div (D(x) ∇v(x)) = 0. ( 34
)
The analogy with diffusion equations makes explicit the direction and magnitude of the smoothing, which correspond respectively to the eigenvectors and eigenvalues of D(x).

The discretization of the gradient and divergence operators yields a large system of equations to be solved. If the system is linear, its sparsity makes it well suited for iterative solvers like Gauss-Seidel or successive over-relaxation (SOR) [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF]. Nevertheless, the linear case is in practice only encountered in the model of [START_REF] Horn | Determining optical flow[END_REF] using quadratic penalties and a linearized data constraint. To cope with non-linearity, the typical approach [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Weickert | Fast parallel algorithms for a broad class of nonlinear variational diffusion approaches[END_REF] is to resort to time-lagged schemes [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] by handling each source of non-linearity with fixed point iterations, turning the problem into a succession of linear systems, and updating iteratively the non-linear parts. Convergence of the scheme is ensured if the linear systems are solved exactly, but approximations and frequent updates often yield in practice good results and much faster convergence.

Fast computational schemes have been employed for achieving near real-time performance. A multigrid framework has been proposed in [START_REF] Bruhn | A confidence measure for variational optic flow methods[END_REF]. The scheme is very efficient, but it is problem-specific and requires a substantial implementation effort. Another approach is to consider the solution of the Euler-Lagrange equation as the steady state of the corresponding diffusion process [START_REF] Bleyer | Segmentation-based motion with occlusions using graph-cut optimization[END_REF], and use the Fast Explicit Diffusion (FED) principle [START_REF] Grewenig | From box filtering to fast explicit diffusion[END_REF] to accelerate convergence by adapting the time steps. The implementation of [START_REF] Gwosdek | A highly efficient GPU implementation for variational optic flow based on the Euler-Lagrange framework[END_REF] exploits the natural parallelization of explicit schemes to achieve a quasi real-time version of the variational method of [START_REF] Zimmer | Optic flow in harmony[END_REF] on GPU, based on FED.

The approach of [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF] has become standard because of its simplicity, the wide range of models it can handle, and its good experimental performances even for non-convex energies [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF][START_REF] Li | Optical flow estimation using Laplacian mesh energy[END_REF][START_REF] Volz | Modeling temporal coherence for optical flow[END_REF].

Early discretization. The computation of the Euler-Lagrange equations can be complex or impossible for some forms of the energy [START_REF] Becker | Variational adaptive correlation method for flow estimation[END_REF]. Moreover, one can argue that the discretization of the Euler-Lagrange equations can generate numerical errors with respect to the original energy [START_REF] Pock | Algorithmic differentiation: Application to variational problems in computer vision[END_REF]. A way to alleviate these shortcomings is to avoid computing the Euler-Lagrange equations and instead to directly discretize the energy [START_REF] Becker | Variational adaptive correlation method for flow estimation[END_REF]. The equation to solve is then simply the cancellation of the differentiated discretized energy:

∂E(w) ∂w = 0. (35)
Solving [START_REF] Boros | Network flows and minimization of quadratic pseudo-boolean functions[END_REF] amounts to inverting a large linear system, similar to the one obtained by Euler-Lagrange equation discretization for energies of the form [START_REF] Black | Learning parameterized models of image motion[END_REF]. Employing a fixed-point iterations scheme to cope with non-linearity amounts to an IRLS approach, which is shown by [START_REF] Liu | Beyond pixels: exploring new representations and applications for motion analysis[END_REF] to be equivalent to the above described resolution of the Euler-Lagrange equations with fixed-point. Contrary to the Euler-Lagrange scheme, the regularization is not imposed to be a derivative of w, and non-local regularization terms (30) can be handled. However, such an approach yields a dense linear system, not solvable with standard iterative methods. In [START_REF] Krähenbühl | Efficient nonlocal regularization for optical flow[END_REF], a linear-time method is proposed to compute matrix product by successive Gaussian filtering, allowing to efficiently inverse the dense linear system with a conjugate gradient solver. The method of [START_REF] Sun | Learning optical flow[END_REF] also exploited the ability of handling more general energies to optimize learned data and regularization terms. The general experimental study of [START_REF] Sun | Secrets of optical flow estimation and their principles[END_REF] and the complex method of [START_REF] Sun | Layered image motion with explicit occlusions, temporal consistency, and depth ordering[END_REF] also follows this approach, with Graduated Non Convexity (GNC) to cope with multimodality of the energy [START_REF] Black | The robust estimation of multiple motions: Parametric and piecewise-smooth flow fields[END_REF][START_REF] Blake | Visual reconstruction[END_REF].

Instead of solving [START_REF] Boros | Network flows and minimization of quadratic pseudo-boolean functions[END_REF], a gradient descent method can be applied to minimize the discretized energy. Due to the large scale of the problem, Newton methods requiring the inversion of the Hessian of the energy are not applicable, and only quasi-Newton methods are computationally tractable. A few works have explored this direction, with Truncated Newton [START_REF] Kalmoun | Trust region versus line search for computing the optical flow[END_REF][START_REF] Kalmoun | Line search multilevel optimization as computational methods for dense optical flow[END_REF] or L-BFGS [START_REF] Dérian | Wavelet-based fluid motion estimation[END_REF].

Half-quadratic minimization. Instead of solving directly the energy [START_REF] Becker | Variational adaptive correlation method for flow estimation[END_REF], a number of optimization methods proceed to the addition of an auxiliary variable splitting the original problem into easier sub-problems. The half-quadratic minimization falls in this class. It can be shown that under non-restrictive conditions [START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF][START_REF] Geman | Constrained restoration and the recovery of discontinuities[END_REF], a function φ can be rewritten as the following function Φ introducing the additional variable γ ∈]0, 1]:

φ(z) = inf γ Φ(z, γ) = inf γ γ z 2 + ψ(γ) , (36)
where the function ψ can be explicitly derived from φ. The robust non-convex penalty φ(z) in data and regularization potentials can be replaced by Φ(z, γ) in the data and regularization terms, so that the optimization in z becomes easy since it involves only quadratic terms. Moreover, the minimization w.r.t. γ has a closed-form solution:

γ = arg min γ Φ(z, γ) = φ (z) 2z . (37)
The original minimization problem in z is thus turned into alternate simple optimizations on z and γ. This approach also leads to the IRLS algorithm described in Section 4.2. The introduction of this approach for optical flow estimation coincided with the use of robust penalties for discontinuity preserving regularization [START_REF] Aubert | Computing optical flow via variational techniques[END_REF][START_REF] Black | Learning parameterized models of image motion[END_REF][START_REF] Deriche | Optical-flow estimation while preserving its discontinuities: a variational approach[END_REF][START_REF] Mémin | Dense estimation and object-based segmentation of the optical flow with robust techniques[END_REF] and was more recently exploited in [START_REF] Corpetti | Dense estimation of fluid flows[END_REF][START_REF] Héas | Bayesian inference of models and hyperparameters for robust optical-flow estimation[END_REF].

Proximal splitting. Another successful optimization method based on alternate minimization of simple sub-problems has been proposed by [START_REF] Aujol | Structure-texture image decomposition -modeling, algorithms, and parameter selection[END_REF] and used for optical flow in [START_REF] Zach | A duality based approach for realtime TV-L1 optical flow[END_REF]. The data and regularization terms are splitted and associated to separate variables, which are quadratically coupled by a third term:

E split (w, v) = Ω ρ data (x, I 1 , I 2 , w)dx E data + 1 2ε Ω w(x) -v(x) 2 dx + λ Ω ρ reg (x, v)dx E reg , (38)
where v is an auxiliary variable. The parameter ε sets the intensity of the coupling. If ε is small, (38) tends to the original univariate problem [START_REF] Becker | Variational adaptive correlation method for flow estimation[END_REF]. The minimization w.r.t. each variable is the computation of the proximal operator of E data and E reg :

       arg min w E split (w, v) = prox E data ( v) arg min v E split ( w, v) = prox E reg ( w) (39)
where the proximity operator of a function f is defined by

prox f (z) = arg min u f (u) + 1 2 u -z 2 . ( 40
)
The minimization problems (39) can be viewed as alternating coarse pixel-wise data matching and denoising of the flow field.

A large number of variants and generalizations of this proximal splitting idea exist, among which the alternating direction method of multipliers (ADMM) [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] or the formalization in a primal-dual framework described in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. This approach has initially been designed for L 1 penalty on the data and regularization terms, a.k.a. T V-L 1 model [START_REF] Aujol | Structure-texture image decomposition -modeling, algorithms, and parameter selection[END_REF][START_REF] Zach | A duality based approach for realtime TV-L1 optical flow[END_REF], for the simplicity of the resulting optimization sub-problems. The optimization of the data term with v fixed is efficiently performed by a thresholding scheme [START_REF] Yang | A fast algorithm for edgepreserving variational multichannel image restoration[END_REF], and fixing w yields the Rudin-Osher-Fatemi model [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], optimized with the duality based algorithm of [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. An advantage is that a differentiable approximation of the L 1 norm is not required, as in [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF], and one can solve for the exact T V-L 1 model.

In a general point of view, the independence of the optimization of data and regularization parts allows one to design dedicated minimization schemes in a variety of cases. The restriction is to be able to compute the proximal operators, and the convergence is ensured only for convex energies. For the data part, the thresholding scheme of [START_REF] Zach | A duality based approach for realtime TV-L1 optical flow[END_REF] is applicable in some cases for the L 1 norm, and efficient solutions have been found to handle an additional fundamental matrix constraint [START_REF] Wedel | Duality TV-L1 flow with fundamental matrix prior[END_REF], a truncated L 1 norm of normalized cross correlation [START_REF] Werlberger | Motion estimation with nonlocal total variation regularization[END_REF], or mutual information [START_REF] Panin | Mutual information for multi-modal, discontinuity-preserving image registration[END_REF]. Another advantage of the pixel-wise nature of the data part minimization is that it enables parallel implementation strategies which can dramatically speed up the algorithm and reach real-time [START_REF] Zach | A duality based approach for realtime TV-L1 optical flow[END_REF]. Based on the pixel-wise property, [START_REF] Steinbrucker | Advanced data terms for variational optic flow estimation[END_REF] proposes a discrete exhaustive matching for optimization in w, which opens the usage of arbitrary data terms, only limited by the computational cost of the matching. Patchbased similarity measures have also recently been implemented with the fast PatchMatch algorithm [START_REF] Barnes | Patchmatch: a randomized correspondence algorithm for structural image editing[END_REF] by [START_REF] Heise | PM-Huber: Patchmatch with Huber regularization for stereo matching[END_REF].

Concerning the regularization part, it is possible to minimize non-local regularization terms [START_REF] Drulea | Motion estimation using the correlation transform[END_REF][START_REF] Werlberger | Efficient minimization of the non-local Potts model[END_REF]. From an evaluation viewpoint, the decoupling in the minization process also allows for a fair comparison of the effects of different data and regularization terms [START_REF] Vogel | An evaluation of data costs for optical flow[END_REF].

Combinatorial methods

In a discrete setting, the solution of the minimization of ( 22) is searched in a set of discrete labels L corresponding to a quantization of the continuous motion vector space or the selection of a finite subset of motion vectors. We give a fast overview of basic principles of optimization methods for this combinatorial problem. For a more complete analysis, see the recent review of [START_REF] Wang | Markov random field modeling, inference & learning in computer vision & image understanding: A survey[END_REF]. The spatial discretization of ( 22) yields an energy in the Markov Random Field (MRF) framework:

E D (w) = x∈Ω D         ρ data (x, I 1 , I 2 , w) + y∈N(x) ρ MRF (x, y, w)         (41)
where Ω D is the discrete image domain, N(x) denotes the pixels interacting with x in the model and ρ MRF is the discrete version of ρ reg (x, y, I 1 , I 2 , w) which explicitly takes into account the interaction between two neighboring pixels x and y.

An important advantage of discrete optimization over the continuous approach is that it does not require differentiation of the energy and can thus handle a wider variety of data and regularization terms. On the counterpart, a trade-off has to be found between the accuracy of the motion labelling and the size of the search space L. Indeed, discrete optimization methods are usually severely limited in terms of accuracy and speed by the number of labels, particularly for optical flow where subpixel accuracy is necessary, and where the two-dimensional motion space becomes more rapidly intractable than the onedimensional stereo case for instance. Therefore, the design of the label space L is a crucial component of discrete optimization methods for optical flow.

Among early methods for minimizing (41), simulated annealing [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF] offers proved convergence towards the global minimum based on stochastic relaxation. However, the optimal convergence is guaranteed only under prohibitive computational cost. An approximate solution can be obtained much faster with the Iterated Conditional Modes [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF] or its High Confidence First variant [START_REF] Chou | The theory and practice of Bayesian image labeling[END_REF], operating by iterative local minimizations of the energy, but this local optimum often yields poor results compared to modern methods [START_REF] Szeliski | A comparative study of energy minimization methods for markov random fields with smoothness-based priors[END_REF], especially for nonconvex functions. A multi-scale MRF approach was designed in [START_REF] Heitz | Multiscale minimization of global energy functions in some visual recovery problems[END_REF] where the energy function is minimized over nested subspaces of the original discrete label space. This constrained optimization scheme offered simple implementation and fast convergence towards high quality estimates for non-linear problems such as motion computation.

Graph cut. The work of [START_REF] Boykov | Markov random fields with efficient approximations[END_REF] gave rise to rapidly growing research interest and advances on graph cut approaches for MRFbased energy minimization. The basis of graph cuts is the maxflow/min-cut problem consisting in finding the optimal path between two nodes in a directed graph, solvable by many algorithm in polynomial time [START_REF] Fulkerson | Flows in networks[END_REF][START_REF] Goldberg | A new approach to the maximum-flow problem[END_REF]. It is possible to model an undirected MRF structure as a directed graph by introducing two additional source and sink nodes, and then interpret the min-cut partition as a binary label segmentation of the MRF energy. The global minimum of the binary optimization can be guaranteed for pairwise interactions and submodular functions. In summary, it is possible to find the global optimum for the energy (41) under the following conditions:

submodularity of ρ MRF (x, y, w), pairwise interactions, binary labels.

Research on graph cuts has attempted to overcome these three constraints, based on the original max-flow/min-cut algorithm.

Submodularity of the pairwise term is a required property for convergence of the algorithm. Finding good approximate solution can be achieved with the Quadratic Pseudo Boolean Optimization (QPBO) [START_REF] Boros | Network flows and minimization of quadratic pseudo-boolean functions[END_REF][START_REF] Kolmogorov | Minimizing nonsubmodular functions with graph cuts-a review[END_REF][START_REF] Rother | Optimizing binary MRFs via extended roof duality[END_REF], leading to an optimal but partial labelling. Dealing with higher-order local interactions between pixels has also been addressed in several recent works [START_REF] Fix | A graph cut algorithm for higher-order Markov random fields[END_REF][START_REF] Ishikawa | Higher-order clique reduction in binary graph cut[END_REF][START_REF] Kohli | P 3 & beyond: Move making algorithms for solving higher order functions[END_REF][START_REF] Komodakis | Beyond pairwise energies: Efficient optimization for higher-order MRFs[END_REF].

Submodularity and pairwise interactions are not too restrictive constraints for a large range of computer vision problems and let room for a number of applications. On the contrary, the cardinality of the label space satisfies |L| > 2 in most cases, so that the binary label requirement is a much harder limitation. The extension of binary graph cuts to multi-label have mostly been realized through iterative move-making techniques. The idea is to create at each iteration a binary-labeled space composed of the current solution and a new proposal label. The label set L is thus explored progressively by each new proposal. If each binary minimization performed with the max-flow/mincut method is ensured to be optimal, a decreasing of the energy is guaranteed at each iteration, and the method converges to a local minimum. It can be noticed that the class of move-making methods is not restricted to graph cuts and the moves can be achieved by other techniques such as the variational approach of [START_REF] Trobin | Continuous energy minimization via repeated binary fusion[END_REF], optimizing in the continuous space by relaxing the binary variable.

The elements of a move-making method are the move-space, specifying the set of new labelling proposals at each iteration, and the way the moves themselves are performed. The most popular move-space is the expansion-move, where the proposal labelling is defined as a constant label map. Another common alternative is the αβ-swap move based on a label exchange at pixels having labels α or β. The range-move [START_REF] Veksler | Graph cut based optimization for MRFs with truncated convex priors[END_REF] extends binary proposal choice to a larger range of labels, in the case of ordered label spaces. Pre-computed labellings computed with independent and possibly continuous methods, can also serve as proposals [START_REF] Lempitsky | Fusion moves for Markov random field optimization[END_REF].

Additionally, computational efficiency of graph-cut approaches have been addressed by several works. The most representative ones are [START_REF] Komodakis | Performance vs computational efficiency for optimizing single and dynamic MRFs: Setting the state of the art with primal-dual strategies[END_REF], operating in a Primal-Dual framework and speeding up convergence by minimizing the Primal-Dual gap, and on the other hand, [START_REF] Kohli | Dynamic graph cuts for efficient inference in Markov random fields[END_REF] working on dynamic MRF and exploiting previous iterations as initializations. For more details about existing move-spaces, see [START_REF] Veksler | Efficient graph-based energy minimization methods in computer vision[END_REF][START_REF] Wang | Markov random field modeling, inference & learning in computer vision & image understanding: A survey[END_REF]. Ap-plications for optical flow have increased rapidly in recent years [START_REF] Chen | Large displacement optical flow from nearest neighbor fields[END_REF][START_REF] Cooke | Two applications of graph-cuts to image processing[END_REF][START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF][START_REF] Glocker | Optical flow estimation with uncertainties through dynamic MRFs[END_REF][START_REF] Glocker | Triangleflow: Optical flow with triangulation-based higher-order likelihoods[END_REF][START_REF] Lempitsky | Fusion moves for Markov random field optimization[END_REF][START_REF] Li | Optical flow estimation using Laplacian mesh energy[END_REF][START_REF] Sun | Layered segmentation and optical flow estimation over time[END_REF][START_REF] Xu | Motion detail preserving optical flow estimation[END_REF].

Message passing. Belief propagation is based on the maxproduct algorithm, which is able to find the MAP of a probability distribution expressed as a product of factors, and thus representable in a factor graph (see [START_REF] Kschischang | Factor graphs and the sum-product algorithm[END_REF] for a detailed introduction). Taking the negative logarithmic version of such distribution amounts to considering MRFs of the form (41), which can motivate to rename the algorithm min-sum in this case. In a nutshell, the minimization is done by iteratively updating local messages reflecting influence of local label configurations on the energy. After convergence of the messages, they can be used to define the probability of assigning a given label to a node, and the label with the maximum probability is chosen. The maxproduct has originally been designed for tree structures and is guaranteed to find the global optimum in this case. Nevertheless, it can also be used for MRFs exhibiting cycles (it is referred to as loopy belief propagation in this case [START_REF] Pearl | Probabilistic reasoning in intelligent systems: networks of plausible inference[END_REF]), without convergence guarantees but showing good experimental results in a large number of computer vision problems [START_REF] Kappes | A comparative study of modern inference techniques for discrete energy minimization problems[END_REF][START_REF] Szeliski | Motion estimation with quadtree splines[END_REF]. The Tree Reweighted message passing approach [START_REF] Wainwright | MAP estimation via agreement on trees: message-passing and linear programming[END_REF] deals with similar concepts, but with a particular message passing strategy based on tree representations. The sequential approach of [START_REF] Kolmogorov | Convergent tree-reweighted message passing for energy minimization[END_REF] has proven to yield good results and computational performance compared to other discrete methods in [START_REF] Szeliski | A comparative study of energy minimization methods for markov random fields with smoothness-based priors[END_REF]. It has been applied for optical flow estimation in [START_REF] Grauer-Gray | Hierarchical belief propagation to reduce search space using cuda for stereo and motion estimation[END_REF][START_REF] Hornácek | Highly overparameterized optical flow using PatchMatch belief propagation[END_REF][START_REF] Lee | Optical flow estimation with adaptive convolution kernel prior on discrete framework[END_REF][START_REF] Mozerov | Constrained optical flow estimation as a matching problem[END_REF]. It is also exploited for image matching involving large deformation models to cope with high computational complexity, as described in [START_REF] Shekhovtsov | Efficient MRF deformation model for non-rigid image matching[END_REF] or in the SIFT-Flow algorithm [START_REF] Liu | Sift flow: Dense correspondence across scenes and its applications[END_REF].

Occlusion handling

The data costs described in Section 3 have been defined as measures of some data constancy between two corresponding points. However, they do not take into account the case of pixels without correspondence, namely occlusions. The occlusion issue in motion estimation, stereovision or image correspondence has motivated different types of investigation [START_REF] Ayvaci | Sparse occlusion detection with optical flow[END_REF][START_REF] Heitz | Multimodal estimation of discontinuous optical flow using Markov random fields[END_REF][START_REF] Bleyer | Segmentation-based motion with occlusions using graph-cut optimization[END_REF][START_REF] Jacobson | An online learning approach to occlusion boundary detection[END_REF][START_REF] Kolmogorov | Computing visual correspondences with occlusion using graph cuts[END_REF][START_REF] Leordeanu | Locally affine sparse-todense matching for motion and occlusion estimation[END_REF][START_REF] Stein | Occlusion boundaries from motion: Lowlevel detection and mid-level reasoning[END_REF][START_REF] Williams | Estimating disparity and occlusions in stereo video sequences[END_REF][START_REF] Xiao | Bilateral filtering-based optical flow estimation with occlusion detection[END_REF]. Occluded regions in the current image of the video sequence are defined as the set of pixels which become hidden by occluding (moving) objects in the next image, as illustrated in Figure 3. These pixels have no corresponding points in the next image and motion is not observable at their locations. Points on (static or moving) objects which leave the camera field of view to the next image are part of a similar situation. Hence, the occlusion issue must be carefully addressed to ensure a reliable and accurate computed optical flow. This is particularly the case for large displacements [START_REF] Strecha | A probabilistic approach to large displacement optical flow and occlusion detection[END_REF], as in the MPI Sintel dataset [START_REF] Butler | A naturalistic open source movie for optical flow evaluation[END_REF]. Most methods implicitly deal with occlusion, together with other sources of violation of the data constancy assumption, by using robust norms in the data term, as emphasized in Section 3. To describe explicit occlusion handling approaches, we must further distinguish between occlusion detection, i.e., segmenting the image into occluded and non-occluded regions, and occlusion filling, i.e., recovering the missing flow subfields in occluded regions. Binary occlusion map (occluded regions in white)

Figure 3: Example of occluded regions.

Occlusion detection

Occlusion detection has been mostly investigated as a subsequent operation to motion computation, by thresholding a consistency measure issued from the estimated motion field. Several criteria have been proposed as geometric forwardbackward motion mismatch [START_REF] Ince | Occlusion-aware optical flow estimation[END_REF][START_REF] Mozerov | Constrained optical flow estimation as a matching problem[END_REF], mapping unicity [START_REF] Xu | Motion detail preserving optical flow estimation[END_REF], divergence of the flow field [START_REF] Ballester | A TV-L1 optical flow method with occlusion detection[END_REF] or data constancy violation [START_REF] Xiao | Bilateral filtering-based optical flow estimation with occlusion detection[END_REF]. The main limitation of the sequential approach is that accuracy of occlusion detection is highly dependent on the quality of the motion estimation. Several flow and image criteria have also been combined in a learning framework [START_REF] Humayun | Learning to find occlusion regions[END_REF].

Other approaches estimate the occlusion map jointly with the motion field in an alternate optimization scheme [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 | Highdimension multilabel problems: Convex or nonconvex relaxation?[END_REF][START_REF] Unger | Joint motion estimation and segmentation of complex scenes with label costs and occlusion modeling[END_REF]. In the continuous setting of [START_REF] Ayvaci | Sparse occlusion detection with optical flow[END_REF] involving a sparsity constraint for occlusion detection, the data constancy deviation is balanced by an estimated continuous residual intensity field, from which occluded points are retrieved by thresholding. The usual drawback of alternate optimization schemes is the convergence towards local minima due to the strong coupling between motion and occlusions. This problem is addressed in [START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF] by introducing an occlusion confidence map in the model, guiding occlusion estimation at each iteration independently from the current motion estimate.

Finally, occlusion map can also be retrieved from the decomposition of the motion field in layers [START_REF] Smith | Layered motion segmentation and depth ordering by tracking edges[END_REF][START_REF] Sun | Layered image motion with explicit occlusions, temporal consistency, and depth ordering[END_REF][START_REF] Sun | Layered segmentation and optical flow estimation over time[END_REF]. Indeed, occluded parts of each layers can be naturally derived from their depth ordering and the way they overlap with each other.

Occlusion filling

When occluded regions are known, occlusion filling with motion vectors is conceptually closely related to image inpainting, since it recovers motion in regions where it is by definition not observable. Inpainting methods can be coarsely divided into two classes, diffusion-based methods [START_REF] Bertalmio | Image inpainting[END_REF][START_REF] Chan | Euler's elastica and curvaturebased inpainting[END_REF] and exemplar-based methods [START_REF] Arias | A variational framework for exemplar-based image inpainting[END_REF][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]. A synthesis of these two approaches has been investigated in [START_REF] Bugeau | A comprehensive framework for image inpainting[END_REF] in a variational framework.

Let us denote O ⊂ Ω the set of occluded pixels of I 1 . Classical methods for occlusion filling with motion vectors operate in a global regularization framework by cancelling the data term at occluded points, that is ρ data (x o , w, I 1 , I 2 ) = 0, ∀x o ∈ O. Thus, this class of diffusion-based (or geometry-oriented) occlusion filling relies on the diffusion process of the regularization to propagate motion from non-occluded regions to occluded regions via partial derivative equation (PDE) resolution [START_REF] Ayvaci | Sparse occlusion detection with optical flow[END_REF][START_REF] Papadakis | Highdimension multilabel problems: Convex or nonconvex relaxation?[END_REF][START_REF] Xu | Motion detail preserving optical flow estimation[END_REF].

It is known that diffusion-based inpainting methods perform well in case of thin missing areas or cartoon-like images, but they are usually unable to handle for large missing regions. Rather than simply cancelling the data term, it is then necessary to design valid data costs for occluded pixels when dealing with large occlusion regions. This is achieved in [START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF], following the inpainting analogy and inspired by the superior performances of exemplar-based approach for large missing regions. At occluded pixels x ∈ O, the data cost imposes proximity with the motion vector of a matched non-occluded pixel m(x o ) ∈ Ω\O, assumed to belong to the same object as x o and thus to have similar motion:

ρ data (x o , w, I 1 , I 2 ) = w(x o ) -w(m(x o )) 2 . ( 42
)
Another possible definition of data cost at occluded pixels can be found in [START_REF] Ballester | A TV-L1 optical flow method with occlusion detection[END_REF], under the assumption of temporal constancy of motion. Assuming that a pixel x o ∈ O of I 1 , occluded in the next frame I 2 , is non-occluded in the previous frame I 0 , backward data constancy is imposed at occluded pixels:

ρ data (x o , w, I 1 , I 2 ) = φ(I 0 -I 1 (x o -w(x o ))). (43)

Combining feature matching and optical flow

An increasing challenge for optical flow estimation is to handle very large displacements and deformations, as it is reflected in the recent MPI Sintel benchmarks [START_REF] Butler | A naturalistic open source movie for optical flow evaluation[END_REF], where is it not rare to find displacements of more than 100 pixels. As a consequence, the gap between feature matching and optical flow tends to vanish, and several methods have tried to combine density and accuracy of optical flow with the ability to capture large displacements of feature matching.

Finding correspondences by matching image features can be considered as a local parametric approach since a given neighborhood of the pixel to match is assumed to translate or undergo another parametric transformation towards its correspondence location. The similarity measures can be patch-based distances (see Section 3.1.2), or more complex and sparse feature descriptions often based on histogram of oriented gradients [START_REF] Bay | Surf: Speeded up robust features[END_REF][START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF][START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF][START_REF] Yu | ASIFT: an algorithm for fully affine invariant comparison[END_REF], or segment matching [START_REF] Wang | MSLD: A robust descriptor for line matching[END_REF]. The fundamental difference lies in the optimization process, since the parametric formulation has a differential optimization process imposing linearization, whereas feature matching explores a discrete space of admissible correspondences. Although being integer displacements and prone to errors, feature matching can thus handle large displacements without coarse-to-fineschemes, with arbitrarily complex similarity measures. Regarding their complementarity of advantages, the combination of feature matching with regularization approaches is therefore of upmost interest.

We consider three approaches to obtain dense and accurate motion fields with feature matching: local filtering of correspondences, integration in a variational framework, and generation of coarse initialization for variational refinement.

Correspondence filtering

Pure feature matching has long been considered unable to produce dense flow fields with competitive accuracy with the previously described local and global approaches. There are three reasons for this:

1. the optimization of the similarity measure is often performed with exhaustive search, which induces prohibitive computational cost, 2. repetitive textures or uninformative regions are sources of ambiguities for the matching process and generate large errors, 3. the correspondence process usually limits the accuracy to integer displacements, contrary to local and global approaches working in the continuous R 2 space.

Several recent advances attempted to overcome these three seemingly inherent limitations.

Research on speeding up block matching include multiscale search strategies [START_REF] Tzovaras | Evaluation of multiresolution block matching techniques for motion and disparity estimation[END_REF], integral images [START_REF] Facciolo | Integral images for block matching[END_REF] or search in trees [START_REF] Kumar | What is a good nearest neighbors algorithm for finding similar patches in images[END_REF], but the recent most spectacular contribution was achieved in [START_REF] Barnes | Patchmatch: a randomized correspondence algorithm for structural image editing[END_REF][START_REF] Barnes | The generalized patchmatch correspondence algorithm[END_REF] with the PatchMatch algorithm. The method scans the image in the lexicographic and inverse order, and alternates two simple steps at each pixel: the propagation step minimizes locally the data cost in the space composed by the current pixel and its two predecessors in the scanning order, the second step proposes a small set of new candidates randomly chosen in the neighborhood of the current correspondence. It is easy to extend the matching to more complex transformations than translation, like rotation or scale factor [START_REF] Barnes | The generalized patchmatch correspondence algorithm[END_REF] by increasing the degree of the search space. The method was originally designed for image editing and was applied to several other applications [START_REF] Barnes | The generalized patchmatch correspondence algorithm[END_REF], with impressive results regarding the low computation time. For motion estimation, the interesting property is that the computational cost is not affected by the spatial extension of the search space, so that no trade-off has to be found between speed and displacement range.

The problem of matching ambiguities comes from the lack of discriminative power of the data cost. Without resorting to explicit regularization (Section 5), the coherency induced by simple local filtering of patch correspondences [START_REF] Hosni | Fast costvolume filtering for visual correspondence and beyond[END_REF] has been shown to be sufficient in practice to reduce most ambiguities and provide satisfactory dense results. The filtering is achieved in [START_REF] Hosni | Fast costvolume filtering for visual correspondence and beyond[END_REF] by guided filtering [START_REF] He | Guided image filtering[END_REF] and in [START_REF] Ma | Constant time weighted median filtering for stereo matching and beyond[END_REF] by weighted median.

Subpixel accuracy is usually reached by upscaling image resolution or by locally approximating the matching surface by a polynomial function. The induced computational cost is reduced in [START_REF] Hosni | Fast costvolume filtering for visual correspondence and beyond[END_REF][START_REF] Steinbrucker | Advanced data terms for variational optic flow estimation[END_REF] by GPU implementation, and can also be handled by iterative refinement [START_REF] Lee | Optical flow estimation with adaptive convolution kernel prior on discrete framework[END_REF].

The combination of the three ingredients has led to the development of competitive optical flow estimation methods based on pure feature matching locally filtered [START_REF] Hosni | Fast costvolume filtering for visual correspondence and beyond[END_REF][START_REF] Ma | Constant time weighted median filtering for stereo matching and beyond[END_REF][START_REF] Tao | Simpleflow: A non-iterative, sublinear optical flow algorithm[END_REF]

Feature matching in global regularized model

As discussed in Section 2.4, a major limitation of the global variational framework is the loss of small and fast objects due to the use of a coarse-to-fine estimation scheme. One recently investigated approach to overcome this problem is to integrate an information from feature matching into the variational framework, thus combining advantages of both methods.

The approach of [START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF], inspired by [START_REF] Héas | Layered estimation of atmospheric mesoscale dynamics from satellite imagery[END_REF], made the first step in this way by adding to the classical data potential a new constraint taking into account an off-line computation of sparse feature correspondences. Let us denote w c the displacement field obtained with a possibly sparse feature matching process. The new combined data potential is then ρ + data (x, I 1 , I 2 , w, w c ) = ρ data (x, I 1 , I 2 , w) + βρ corresp (x, w, w c ) (44) where β is a trade-off parameter and the matching potential is defined as: ρ corresp (x, w, w c ) = δ(x, w c ) c(x, w c ) φ( w(x)w c (x) 2 ). [START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF] The binary function δ(x, w c ) returns 1 if w c (x) is defined at x and 0 otherwise, and the weights c(x, w c ) correspond to the matching cost. The third term of (45) imposes the motion field to be close to the motion vectors obtained by feature matching.

The advantage of this approach is that the term ρ corresp (x, w, w c ) is both differentiable and valid for large displacements. Problems related to the use of coarse-to-fine schemes are thus avoided. The main drawback is that the importance of the matching term ρ corresp (x, w c , w) relatively to the data term ρ data (x, I 1 , I 2 , w) is mostly determined by the value of the matching cost c(x, w c ). Consequently the final variational estimation is extremely dependent on the reliability of the confidence measure, which is very difficult to guarantee as emphasized in [START_REF] Braux-Zin | A general dense image matching framework combining direct and feature-based costs[END_REF]. Matching errors are thus easily driven by c(x, w c ) and are likely to have a high impact on the final result. Reducing the impact of local feature matching errors by the regularization and the robust penalization φ(•) is insufficient in practice in a lot of cases.

Recently, the authors of [START_REF] Weinzaepfel | Deepflow: Large displacement optical flow with deep matching[END_REF] based their method on [START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF] by taking the model [START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF] and improving the feature matching stage. They demonstrated a significant performance improvement due to the increased reliability of the matching. The method of [START_REF] Braux-Zin | A general dense image matching framework combining direct and feature-based costs[END_REF] is also built upon [START_REF] Brox | Large displacement optical flow: descriptor matching in variational motion estimation[END_REF] with a modified matching component by using segment matching [START_REF] Wang | MSLD: A robust descriptor for line matching[END_REF]. The matching term is generalized to handle weakly localized line features. It is done through a point-to-line distance in addition to the pointto-point distance of (45), combined with a confidence measure for segment matching based on the assumption of a linear mapping between segments. Moreover, considering the unreliability of confidence measure, [START_REF] Braux-Zin | A general dense image matching framework combining direct and feature-based costs[END_REF] gets rid of the matching cost c(x, w c ) and relies only on the regularization for robustness to matching errors. To sum up, the tendency is to take into account the great dependency of the estimation on the quality of the feature matching, and thus to concentrate most efforts on the design of robust matching algorithms.

Feature matching for coarse initialization

Another recently investigated approach taking advantage of feature matching consists in exploiting the (possibly sparse) field of integer displacements to provide a coarse but relevant initialization for a variational refinement [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]. These methods are composed of three steps:

1. a feature matching step provides a limited number of candidates at each pixel;

2. these candidates serve as labels for the discrete optimization (see Section 5.2) of a global regularized model;

3. the resulting coarse optical flow field is refined with one of the variational optimization techniques described in Section 5.2.

The idea is that steps 1 and 2 handle indifferently small and large displacements of objects at any scale, and the initialization is assumed to be good enough to avoid the coarse-to-fine scheme in step 3. The specificity of [START_REF] Xu | Motion detail preserving optical flow estimation[END_REF] is that the three steps are repeated at each level of a multiresolution pyramid. Since steps 2 and 3 have already been discussed in previous subsections, we focus on the specificity of step 1, that is, the production of candidates from feature matching.

In [START_REF] Xu | Motion detail preserving optical flow estimation[END_REF], the feature matching is only performed at a restricted number of keypoints. After pruning of similar vectors, only a few displacement vectors are retained. Each of these vectors is expanded to produce a global constant flow field, used as candidates for step 2.

The approach of [START_REF] Mozerov | Constrained optical flow estimation as a matching problem[END_REF] is to consider an integer discretization of the two-dimensional motion field. Based on the observation that correlation-based patch matching is able to reproduce coarsely the motion distribution pattern of ground truth motion fields, the discretized motion space is delineated by the matching vectors.

The approach of [START_REF] Chen | Large displacement optical flow from nearest neighbor fields[END_REF] is close to [START_REF] Mozerov | Constrained optical flow estimation as a matching problem[END_REF] since the idea is to rely on the dominant motion patterns of dense patch matching. In [START_REF] Chen | Large displacement optical flow from nearest neighbor fields[END_REF], the candidates to feed discrete optimization are obtained by explicitly clustering PatchMatch motion fields [START_REF] Barnes | The generalized patchmatch correspondence algorithm[END_REF] to keep only dominant patterns.

Another strategy has recently been investigated in [START_REF] Fortun | Aggregation of local parametric candidates with exemplar-based occlusion handling for optical flow[END_REF], exploiting the aggregation paradigm. Feature matching is used as a part of a generation procedure of continuous motion candidates, and the best candidates are then selected in the aggregation step with a global model, without need of subsequent variational refinement.

Conclusion and perspectives

Through the analysis of optical flow literature, we have introduced principles of optical flow modeling and computation, and classified the main methodological aspects of existing methods. From this taxonomy, interests and limitations of each class have been put forward.

This study also exhibits the main research directions for improving current state-of-the-art. In terms of performance, it can be considered from results of the Middlebury benchmark [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF] that for a range of sequences with small displacements, moderate intensity changes and piecewise smooth displacements of large structures, globally regularized models and joint estimation and segmentation methods are able to obtain satisfying results. However, moving away from these dataset characteristics, as it is the case in recent benchmarks [START_REF] Butler | A naturalistic open source movie for optical flow evaluation[END_REF][START_REF] Geiger | Are we ready for autonomous driving? the KITTI vision benchmark suite[END_REF] or applicative domains [START_REF] Fortun | Aggregation of patch-based estimations for illumination-invariant optical flow in live cell imaging[END_REF], unveils the shortcomings of current methods in challenging situations, and it opens the way to potentially large improvements.

In particular, some of the main difficulties are caused by large displacements. To avoid the drawbacks of coarse-to-fine schemes of variational methods, an active research direction focuses on the design of computationally tractable discrete optimization methods for the large scale problem of optical flow. We have noticed that dealing with larger displacements brings optical flow closer to feature matching. However, attempts to combine these two worlds still consider them separately and try to fuse their independent results, often leading to high sensitivity to matching errors. Another essential problem lies in large occlusions that are inevitably associated to large motion. Large intensity changes appearing in real environments or in specific applications, or caused by large deformations, is another main problem of current methods. The most promising direction seems to be the spatial adaptation of the data term.

Finally, another remaining challenge is to reduce the computational cost of the best performing methods.

We observe an evolution of the scope of optical flow, which was originally devoted to image sequences with high frame-rate, for which deviations from small displacements and brightness conservation stayed reasonably low. As this well identified problem becomes better handled, optical flow estimation tends to apply on larger class of sequences. Thus, the gap between optical flow and other correspondence problems like registration or feature matching progressively reduces and opens the way to unified correspondence methods.

Figure 1 :

 1 Figure 1: Two types of visualization of the motion field transforming I 1 in I 2 .

Figure 2 :

 2 Figure 2: Examples of frames and groud truth from the main optical flow evaluation benchmarks.

  investigated higher-order data constancy like Laplacian of the image f (I) = ∆I, or norm and determinant of the Hessian f (I) = |H I|, f (I) = det(H I).

I 1 I 2

 12 Occlusions on I 1 (occluded regions in red)

http://www.cvlibs.net/datasets/kitti/index.php

Acknowledgements

This work was realized as part of the Quaero program, funded by OSEO, French State agency for innovation. It was also partly supported by the France-BioImaging project granted by the Investissement d'Avenir program, and by UMR 144 CNRS of Institut Curie.