Robust optical flow estimation

Published in Image Processing On Line on 2013–10–28. Submitted on 2012–06–22, accepted on 2013–05–31. ISSN 2105–1232 c 2013 IPOL & the authors CC–BY–NC–SA This article is available online with supplementary materials, software, datasets and online demo at http://dx.doi.org/10.5201/ipol.2013.21 2013/07/19 v0.4 IPOL article class Robust Optical Flow Estimation Javier S´anchez1, Nelson Monz´on2and Agust´ın Salgado3 1CTIM, University of Las Palmas de Gran Canaria, Spain ([email protected]) 2CTIM, University of Las Palmas de Gran Canaria, Spain ([email protected]) 3CTIM, University of Las Palmas de Gran Canaria, Spain ([email protected]) Abstract In this work, we describe an implementation of the variational method proposed by Brox et al. in 2004, which yields accurate optical flows with low running times. It has several benefits with respect to the method of Horn and Schunck: it is more robust to the presence of outliers, produces piecewise-smooth flow fields and can cope with constant brightness changes. This method relies on the brightness and gradient constancy assumptions, using the information of the image intensities and the image gradients to find correspondences. It also generalizes the use of continuous L1functionals, which help mitigate the effect of outliers and create a Total Variation (TV) regularization. Additionally, it introduces a simple temporal regularization scheme that enforces a continuous temporal coherence of the flow fields. Source Code The source code, the code documentation, and the online demo are accessible at the IPOL web page of this article1. In this page an implementation is available for download. This file contains two directories: one for the spatial method and another for the temporal method. The spatial method is suitable for general image sequences, while the temporal method should be used when the flow fields are known to be very continuous. Keywords: optical flow, motion estimation, variational techniques, PDE 1 Introduction The estimation of accurate motion fields is an important challenge in image processing and computer vision. Among the most accurate methods in the literature, variational approaches have proven an outstanding performance with respect to other strategies. Starting from the first variational model of Horn and Schunck [5], we find an extensive literature on such methods. Brox et al. [3] proposed a technique that is based on differentiable L1functionals. This kind of functionals allows to create piecewise-smooth flow fields. A more efficient TV-L1optical flow method was presented in Zach et al. [11], which has also been implemented on an IPOL article [9]. The 1http://dx.doi.org/10.5201/ipol.2013.21 Javier S´anchez, Nelson Monz´on and Agust´ın Salgado, “Robust Optical Flow Estimation,” Image Processing On Line, vol. 2013, pp. 242–260. http://dx.doi.org/10.5201/ipol.2013.21 Robust Optical Flow Estimation benefit of the Brox et al. method is that it keeps a better spatial coherence between both unknowns of the optical flow, at the expense of creating rounded effects at flow discontinuities. The Zach et al. method creates sharper discontinuities, but it strongly suffers from the staircase effect of pure TV schemes. Rounded effects at motion discontinuities can be mitigated with the use of decreasing functions, like in the work of Wedel et al. [10]. A different alternative is to use image-based anisotropic schemes, like in the work of ´ Alvarez et al. [1], that allows respecting the image discontinuities. In this work, we implement the spatial and temporal methods as in the original proposal. Both methods are very similar, with the main difference being that the temporal method includes a continuous smoothing scheme in the temporal dimension. This is suitable when the optical flow functions are smooth. Nevertheless, the flow discontinuities are typically degraded by the temporal regularizer. In the presence of large displacements, it is better to use a nonlinear temporal scheme, like in Salgado and S´anchez [8]. Given a sequence of images, I: Ω ⊂R3→R,of gray level values in space and time, x= (x, y, t)T∈Ω, the optical flow is defined as a dense mapping, w= (u(x), v(x),1)T, between the pixels of every two consecutive images. The scalar fields u(x) and v(x) are the xand ydisplacements in the 3D volume, respectively. Each frame is considered to be at a distance 1 in time from the previous and following frames. The spatial gradient of the image is given by ∇I= (Ix, Iy)T, with Ix, Iythe first order derivatives in xand y. The gradient of the optical flow is defined as ∇u=((ux, uy)T, (ux, uy, ut)T.(1) The first spatial gradient is used in the Spatial method, where the optical flow is computed between two frames. The second spatio-temporal gradient is used in the Temporal method. Since the abstract framework is identical in both cases, we keep the same notation for both methods. This distinction appears later in the numerical scheme (section 2). We suppose that the pixel intensities remain constant along the trajectories of the moving particles. This is normally referred to as the brightness constancy assumption, which in a continuous setting yields the optical flow constraint equation, proposed by Horn and Schunck [5] (see also the related IPOL implementation [6]). In general, we suppose that our problem is not continuous, so that the brightness constancy assumption can be written in its nonlinear form as I(x+w)−I(x)=0.(2) The following energy model is a slight variation of the one proposed by Brox et al. in 2004 [3]: E(w) = Z Ω Ψ(I(x+w)−I(x))2dx +γZ Ω Ψ|∇I(x+w)−∇I(x)|2dx +αZ Ω Ψ|∇u|2+|∇v|2dx.(3) with Ψ(s2) = √s2+2and := 0.001 a small constant. In the original proposal [3], the brightness and gradient constancy terms (first two terms) were included inside the same Ψ function. As proposed by Bruhn and Weickert [4], the separation of these two assumptions, like in equation (3), is better justified. In this sense, our model is more similar to the one by Bruhn and Weickert [4]. The last term corresponds to the regularizing strategy, which is a differentiable approximation of the TV scheme. The method of Zach et al. [11], or Wedel et al. [10], uses a pure TV scheme and 243 Javier S´ anchez, Nelson Monz´ on and Agust´ ın Salgado separates both components of the optical flow. Although both approaches rely on L1functionals, they provide very different results, as can be seen in the IPOL article of S´anchez et al. [9]. The minimum of the previous energy model can be found by solving the associated Euler-Lagrange equations, given by 0 =Ψ0 D·(I(x+w)−I(x)) ·Ix(x+w) +γΨ0 G·((Ix(x+w)−Ix(x)) ·Ixx(x+w)+(Iy(x+w)−Iy(x)) ·Ixy(x+w)) −αdiv (Ψ0 S·∇u), 0 =Ψ0 D·(I(x+w)−I(x)) ·Iy(x+w) +γΨ0 G·((Ix(x+w)−Ix(x)) ·Ixy(x+w)+(Iy(x+w)−Iy(x)) ·Iyy(x+w)) −αdiv(Ψ0 S·∇v),(4) with Ψ0(s2) = 1 2√s2+2. In order to simplify the equations, we use the following notation: Ψ0 D:=Ψ0(I(x+w)−I(x))2, Ψ0 G:=Ψ0|∇I(x+w)−∇I(x)|2, Ψ0 S:=Ψ0|∇u|2+|∇v|2.(5) Equation (4) cannot be solved easily because it is nonlinear in expressions like I(x+w). To avoid these nonlinearities, we use first order Taylor expansions. We introduce an index, k, so that our current unknown, wk+1, depends on a previous constant value, wk. I(x+wk+1)≈I(x+wk) + Ix(x+wk)duk+Iy(x+wk)dvk Ix(x+wk+1)≈Ix(x+wk) + Ixx(x+wk)duk+Ixy(x+wk)dvk Iy(x+wk+1)≈Iy(x+wk) + Ixy(x+wk)duk+Iyy(x+wk)dvk,(6) with wk= (uk, vk)T,duk=uk+1 −ukand dvk=vk+1 −vk. We assume that wkis a close approximation to our unknown wk+1. The original method, and our implementation, works directly with the motion increments (duk, dvk). The optical flow is incrementally updated from the motion increment as uk+1 =uk+duk,vk+1 =vk+dvk. This is the strategy followed in other works like, for instance, in M´emin and P´erez [7]. A different alternative is to compute the full optical flow directly, like in the work by ´ Alvarez et al. [1]. There still remains another nonlinearity due to the Ψ0functions. Thus, our numerical scheme should be enclosed in two fixed point iterations: the outer iterations, k, related with the stability of the Taylor expansions; and the inner iterations, l, that account for the nonlinearities of the Ψ0 functions. Therefore, the system of equations reads as 0 =(Ψ0 D)k,l ·I(y) + Ix(y)duk,l+1 +Iy(y)dvk,l+1 −I(x)·Ix(y) +γ(Ψ0 G)k,l ·Ix(y) + Ixx(y)duk,l+1 +Ixy(y)dvk,l+1 −Ix(x)·Ixx(y) +Iy(y) + Ixy(y)duk,l+1 +Iyy(y)dvk,l+1 −Iy(x)·Ixy(yk,l) −αdiv (Ψ0 S)k,l ·∇(uk,l +duk,l+1) 0 =(Ψ0 D)k,l ·I(y) + Ix(y)duk,l+1 +Iy(y)dvk,l+1 −I(x)·Ix(y) +γ(Ψ0 G)k,l ·Ix(y) + Ixx(y)duk,l+1 +Ixy(y)dvk,l+1 −Ix(x)·Ixx(y) +Iy(y) + Ixy(y)duk,l+1 +Iyy(y)dvk,l+1 −Iy(x)·Ixy(yk,l) −αdiv (Ψ0 S)k,l ·∇(vk,l +dvk,l+1),(7) with y=x+wk,l. 244 Robust Optical Flow Estimation 2 Numerical Scheme The system of equations (7) can be efficiently solved using the SOR method. The unknowns duk,l+1 and dvk,l+1, in pixel (i, j, k), are expressed in terms of the remaining terms, and their values are iteratively updated until the method converges to a steady state solution. In this sense, we introduce an additional fixed point iteration scheme, s, for the SOR method. Partial derivatives are approximated using central differences. The discretization of the divergence is separated in three variables: div (Ψ0 S)k,l ·∇(uk,l +duk,l+1)= div (Ψ0 S)k,l ·∇uk,l+ div (Ψ0 S)k,l ·∇duk,l+1≈div u + (div du −div d ·duk,l+1 i,j,k ), where div u discretizes the first divergence term, div d and div du correspond to the second term. In the second term, div du stands for the values corresponding to the neighbors of du, and div d stands for the coefficients accompanying du at the current pixel, duk,l+1 i,j,k . These variables are given by the following expressions: div u :=(Ψ0 S)i+1,j,k + (Ψ0 S)k,l i,j,k 2uk,l i+1,j,k −uk,l i,j,k+(Ψ0 S)i−1,j,k + (Ψ0 S)k,l i,j,k 2uk,l i−1,j,k −uk,l i,j,k+ (Ψ0 S)i,j+1,k + (Ψ0 S)k,l i,j,k 2uk,l i,j+1,k −uk,l i,j,k+(Ψ0 S)i,j−1,k + (Ψ0 S)k,l i,j,k 2uk,l i,j−1,k −uk,l i,j,k+ (Ψ0 S)i,j,k+1+ (Ψ0 S)k,l i,j,k 2uk,l i,j,k+1−uk,l i,j,k+(Ψ0 S)i,j,k−1+ (Ψ0 S)k,l i,j,k 2uk,l i,j,k−1−uk,l i,j,k.(8) div du :=(Ψ0 S)i+1,j,k + (Ψ0 S)k,l i,j,k 2duk,l+1 i+1,j,k +(Ψ0 S)i−1,j,k + (Ψ0 S)k,l i,j,k 2duk,l+1 i−1,j,k+ (Ψ0 S)i,j+1,k + (Ψ0 S)k,l i,j,k 2duk,l+1 i,j+1,k +(Ψ0 S)i,j−1,k + (Ψ0 S)k,l i,j,k 2duk,l+1 i,j−1,k+ (Ψ0 S)i,j,k+1+ (Ψ0 S)k,l i,j,k 2duk,l+1 i,j,k+1+(Ψ0 S)i,j,k−1+ (Ψ0 S)k,l i,j,k 2duk,l+1 i,j,k−1.(9) div d :=(Ψ0 S)i+1,j,k + (Ψ0 S)k,l i,j,k 2+(Ψ0 S)i−1,j,k + (Ψ0 S)k,l i,j,k 2+ (Ψ0 S)i,j+1,k + (Ψ0 S)k,l i,j,k 2+(Ψ0 S)i,j−1,k + (Ψ0 S)k,l i,j,k 2+ (Ψ0 S)i,j,k+1+ (Ψ0 S)k,l i,j,k 2+(Ψ0 S)i,j,k−1+ (Ψ0 S)k,l i,j,k 2.(10) These expressions are the same for the other component of the optical flow, changing uby v. The last two terms of equations (8), (9) and (10), in italics, correspond to the temporal regularization of the optical flow. This is implemented in the Temporal method and removed in the Spatial one. The finite difference scheme in the temporal method is computed using information from the previous, k−1, and following, k+ 1, frames. If we define y=x+wk,l and separate the parts of the equation 245 Javier S´ anchez, Nelson Monz´ on and Agust´ ın Salgado that remain constant during the SOR iterations, we may define the following variables: Au := −(Ψ0 D)k,l (I(y)−I(x)) Ix(y) + α div u, −γ(Ψ0 G)k,l ((Ix(y)−Ix(x)) Ixx(y)+(Iy(y)−Iy(x)) Ixy(y)) , Av := −(Ψ0 D)k,l (I(y)−I(x)) Iy(y) + α div v −γ(Ψ0 G)k,l ((Ix(y)−Ix(x)) Ixy(y)+(Iy(y)−Iy(x)) Iyy(y)) , Du :=(Ψ0 D)k,lI2 x(y) + γ(Ψ0 G)k,l I2 xx(y) + I2 xy(y)+α div d, Dv :=(Ψ0 D)k,lI2 y(y) + γ(Ψ0 G)k,l I2 yy(y) + I2 xy(y)+α div d, D:=(Ψ0 D)k,lIx(y)Iy(y) + γ(Ψ0 G)k,l (Ixx(y) + Iyy(y)) Ixy(y).(11) In order to compute expressions like I(x+wk,l), we use bicubic interpolation. Putting all together, we arrive to the SOR scheme, which is given by duk,l,s+1 :=(1 −w)duk,l,s +wAu −D·dvk,l,s+1 +α div du Du , dvk,l,s+1 :=(1 −w)dvk,l,s +wAv −D·duk,l,s+1 +α div dv Dv ,(12) with w∈(0,2) the SOR relaxation parameter. In our implementation, we choose w= 1.9 by default. This numerical approximation is calculated until the method converges to a steady state solution or it exceeds a maximum number of iterations. The stopping criterion is 1 NX i,j,k dus+1 i.j,k −dus i.j,k2+dvs+1 i.j,k −dvs i.j,k2< ε2,(13) with Nthe number of pixels in all frames and εthe stopping criterion threshold. The iterative process stops when condition (13) is true or a maximum number of iterations is reached. This is different from the original article [3], where a fixed number of SOR iterations is used. Once it has converged, we go to the next inner iteration, l+ 1, and restart the variables in (11). In our implementation, we make use of the OpenMP library to enable multiple threads during the computations. This allows us to take advantage of multiple cores, accelerating the run time of the algorithm. Nevertheless, due to the nature of the SOR scheme, it is probable that the results are slightly different, even when the process is launched with the same number of processors. 3 Pyramidal Structure In order to estimate large displacements, we embed the optical flow method in a pyramidal structure. We follow the same strategy presented in a previous IPOL article [6] and reproduce here the basic ideas. Our algorithm creates a pyramid of down-sampled images. The pyramid is created by reducing the images by a factor η∈(0,1). Before downsampling, the images are smoothed with a Gaussian kernel of a standard deviation that depends on η. For a set of scales s= 0,1, . . . , Nscales −1, the pyramid of images is built as Is(ηx) := Gσ∗Is−1(x).(14) After the convolution, the images are sampled using bicubic interpolation. The value of σdepends on ηand is calculated as σ(η) := σ0pη−2−1, with σ0:= 0.6.(15) 246 Robust Optical Flow Estimation Then, starting at the coarsest scale, the system of equations is solved in each scale to get successive approximations of the optical flow. Every intermediate solution is used as initialization in the following scale. To transfer the values from a coarser scale, the flow field is updated as us−1(x) := 1 ηus(ηx) vs−1(x) := 1 ηvs(ηx) (16) 4 Parameters of the Method This method depends on the parameters given in table 1. These parameters are: αand γ, which define the smoothness and the preservation of gradient structures in the optical flow, respectively; the parameters for the pyramidal scheme, Nscales and η, that stand for the number of scales and the downsampling factor; and the parameters for the numerical scheme, composed of the inner and outer iterations, and the stopping criterion threshold (ε). Table 1: Parameters of the method Parameter Explanation αRegularization parameter. It determines the smoothness of the output. The bigger this parameter is, the smoother the solutions we obtain. γParameter associated with the gradient constancy term in equation (3). Nscales Number of scales in the pyramidal structure. If the flow field is very small (about one pixel), it can be set to 1. Otherwise, it should be set so that (1/η)N−1is larger than the expected size of the largest displacement (see the previous article by Meinhardt-Llopis and S´anchez [6] for more details). ηDownsampling factor. It is used to downscale the original images in order to create the pyramidal structure. Its value must be in the interval (0,1). With η= 0.5, the images are reduced to half their size in each dimension from one scale to the following. εStopping criterion threshold. It is the threshold used to stop the SOR iterations, given in equation (13). inner iterations Number of inner iterations in the numerical scheme. It corresponds to index lin equation (12). outer iterations Number of outer iterations in the numerical scheme. It corresponds to index kin equation (12). 5 Algorithm Next, we describe the algorithm that implements the numerical scheme in equation (12). The algorithm takes a set of gray level images as input data and computes the optical flows between every pair of consecutive images. The number of calculated optical flows is one less than the number of input images. We separate the algorithm in two modules: one procedure that computes the optical flows in each scale, and the main algorithm that is in charge of handling the pyramidal structure. In the procedure, MAXITER is the maximum number of iterations allowed for the convergence of the SOR method. Its value is constant and is high enough to let the method converge. In the 247 Javier S´ anchez, Nelson Monz´ on and Agust´ ın Salgado Procedure brox optic flow(I, u, v, α, γ, ε, inner iterations, outer iterations) Compute Ix, Iy 1 Compute Ixx, Iyy, Ixy 2 for no ←0to outer iterations −1do3 Compute I(x+w), Ix(x+w), Iy(x+w) using bicubic interpolation4 Compute Ixx(x+w), Ixy(x+w), Iyy(x+w) using bicubic interpolation5 Compute ux, uy, vx, vy 6 Compute Ψ0 Susing equation (5)7 Compute div u, div v, div d using equations (8) and (10)8 du ←09 dv ←010 for ni ←0to inner iterations −1do11 Compute Ψ0 D,Ψ0 Gusing equation (5)12 Compute Au, Av, Du, Dv, D using equation (11)13 while error > ε and nsor < MAXITER do14 du ←(1 −ω)du +ω(Au −D dv +α div du)/Du15 dv ←(1 −ω)dv +ω(Av −D du +α div dv)/Dv16 Compute error with equation (13)17 nsor ←nsor + 118 end19 end20 u←u+du21 v←v+dv22 end23 Algorithm 2: Pyramidal structure management Input:I, u, v, α, γ, Nscales, η, ε, inner iterations, outer iterations Output:u, v Normalize images between 0 and 2551 Convolve the images with a Gaussian of σ= 0.82 Create the pyramid of images Isusing η(with s= 0, . . . , Nscales −1)3 for s←Nscales −1to 0do4 brox optic flow(I, us, vs, α, γ, inner iterations, outer iterations)5 if s > 0then6 us−1(x) := 1 ηus(ηx) 7 vs−1(x) := 1 ηvs(ηx) 8 end9 end10 248 Robust Optical Flow Estimation source code, the SOR loop is unrolled in order to avoid boundary tests when computing div(duk+1) and div(dvk+1). This means that the first and last columns and rows, and the four corners of the images, are computed separately. The procedure for the spatial and temporal methods is basically the same, except for the computation of the divergence terms. They differ in the calculation of variables Ψ0 S,div u,div v and div d. In theory, as we have seen before, the method needs three levels of iterations. Nevertheless, we have found in the experiments that the inner iterations can be integrated in the outer iterations without a loss of precision. For this reason we have decided to estimate the values of Ψ0 S,div u, div v and div d inside the outer iterations. We prefer this option in order to avoid the estimation of ∇du. The other computations, at the beginning of the inner loop, may also be integrated in the outer loop. The main process is given in algorithm 2. In order to turn the method more stable to the input parameters, it first normalizes the images between 0 and 255; it convolves the finest scale images with a small Gaussian kernel; then, it creates the pyramidal structure for the whole sequence; and, finally, it goes over the set of scales computing the optical flows at different resolutions. 6 Experimental Analysis In this section, we examine the behavior of the method for some standard image sequences. We have used the Yosemite and Yosemite with Clouds sequences, which have also been analyzed by Brox et al. in 2004 [3]. On the other hand, we use the RubberWhale and Urban2 sequences from the Middlebury benchmark database2[2]. The results are shown in figure 2(the color scheme used to represent the orientation and magnitude of optical flows is displayed in figure 1). In these experiments, αand γare adapted to get the best results for each sequence. The values for the remaining parameters are: η= 0.75, ε= 0.0001, inner iterations = 1, outer iterations = 38, and Nscales is automatically calculated so that the coarsest scale works with images around 16 ×16 pixels. Figure 1: Color scheme used to represent the orientation and magnitude of optical flows. In general, the method respects the motion discontinuities and creates piecewise continuous flow fields. This is due to the L1functionals and the TV regularization scheme. In the Yosemite with Clouds sequence, the sky motion is translational (2 pixels to the right). Although there are illumination changes in the clouds, the solution obtained by the method is very accurate. This is due to the gradient constancy term. Note that γis larger than in Yosemite. This shows that the method can correctly handle constant brightness shifts, provided that the gradient does not vary. Tables 2and 3show the Average Angular Error (AAE) and Average End-point Error (EPE) for the Spatial and Temporal methods, respectively. The AAE and EPE are calculated as in the 2http://vision.middlebury.edu/flow/ 249 Javier S´ anchez, Nelson Monz´ on and Agust´ ın Salgado Sequence Ground truth Spatial Temporal Yosemite Yosemite with Clouds RubberWhale Urban2 Figure 2: Results for the Yosemite, Yosemite with Clouds, RubberWhale and Urban2 sequences. First column shows frame 6 for Yosemite and Yosemite with Clouds, and frame 10 for RubberWhale and Urban2. Second column shows the corresponding ground truth optical flows. Third and fourth columns show the results for the Spatial and Temporal methods, respectively. α= 2 α= 7 α= 15 α= 30 Figure 3: Results for the Urban2 sequence using the temporal method. 250 Robust Optical Flow Estimation 7 Examples In this section, we show the results for the sequences in the Middlebury benchmark database [2]. Figure 12 depicts the results for all the tests sequences in the database, except Dimetrodon and Venus for the temporal method. This is because there are not enough images to use the method properly. Table 5show the AAE and EPE obtained for α= 18 and γ= 7 in the spatial method, and α= 2.5, γ= 2 in the temporal method. In both cases, we have set the following parameters: η= 0.75, ε= 0.0001, inner iterations = 1 and outer iterations = 15. Nscales is automatically calculated so that the coarsest scale works with images around 16 x 16 pixels. Table 5: AAE and EPE for the Middlebury test sequences. Error Grove2 Grove3 Hydrangea RubberWhale Urban2 Urban3 Dimetrodon Venus Spat. AAE 2.455o6.481o2.442o3.696o2.561o4.804o1.663o4.599o Spat. EPE 0.174 0.693 0.200 0.111 0.368 0.544 0.086 0.292 Temp. AAE 2.569o7.031o4.468o5.435o5.90o3 6.681o- - Temp. EPE 0.184 0.796 0.346 0.168 0.628 0.784 - - Finally, in figure 13 we show the results for the evaluation sequences using the same parameter configuration (spatial method). 8 Video This is an example of applying the spatial optical flow method, frame by frame, to a video: KarlWilhelm-Straβe3(21,1 MB). The original video of this traffic sequence can be found at the Institut f¨ur Algorithmen und Kognitive Systeme4web pages. Acknowledgements This work has been partially supported by the Spanish Ministry of Science and Innovation through the research project TIN2011-254885. Image Credits All images by the authors except: Daniel Scharstein, Middlebury benchmark database6. Lynn Quam. 3http://dx.doi.org/10.5201/ipol.2013.21/data/dt_passat03.avi 4http://i21www.ira.uka.de/image_sequences/ 5http://www.ctim.es/projects/opticalflow2012/ 6http://vision.middlebury.edu/flow/data/ 257 Javier S´ anchez, Nelson Monz´ on and Agust´ ın Salgado Sequence Ground truth Spatial method Temporal method Figure 12: Results for the Middlebury test sequences. 258 Robust Optical Flow Estimation Sequence Optical flow Sequence Optical flow Figure 13: Results for the Middlebury evaluation sequences. 259 Javier S´ anchez, Nelson Monz´ on and Agust´ ın Salgado References [1] Luis ´ Alvarez, Joachim Weickert, and Javier S´anchez. Reliable estimation of dense optical flow fields with large displacements. International Journal of Computer Vision, 39(1):41–56, 2000. http://dx.doi.org/10.1023/A:1008170101536. [2] Simon Baker, Daniel Scharstein, J. P. Lewis, Stefan Roth, Michael J. Black, and Richard Szeliski. A database and evaluation methodology for optical flow. In International Conference on Computer Vision, pages 1–8, 2007. http://dx.doi.org/10.1109/ICCV.2007.4408903. [3] Thomas Brox, Andr´es Bruhn, Nils Papenberg, and Joachim Weickert. High accuracy optical flow estimation based on a theory for warping. In T. Pajdla and J. Matas, editors, European Conference on Computer Vision (ECCV), volume 3024 of Lecture Notes in Computer Science, pages 25–36, Prague, Czech Republic, May 2004. Springer. http://dx.doi.org/10.1007/ 978-3-540-24673-2_3. [4] Andr´es Bruhn and Joachim Weickert. Towards ultimate motion estimation: Combining highest accuracy with real-time performance. In International Conference on Computer Vision (ICCV), volume 1, pages 749–755, Washington, DC, USA, October 2005. IEEE Computer Society. http: //dx.doi.org/10.1109/ICCV.2005.240. [5] Berthold K. P. Horn and Brian G. Schunck. Determining optical flow. Artificial Intelligence, 17:185–203, 1981. http://dx.doi.org/10.1016/0004-3702(81)90024-2. [6] Enric Meinhardt-Llopis, Javier S´anchez P´erez, and Daniel Kondermann. Horn-Schunck Optical Flow with a Multi-Scale Strategy. Image Processing On Line, 2013:151–172, 2013. http: //dx.doi.org/10.5201/ipol.2013.20. [7] E. M´emin and P. P´erez. Dense estimation and object-based segmentation of the optical-flow with robust techniques. IEEE Transactions on Image Processing, 7(5):703–719, May 1998. http://dx.doi.org/10.1109/83.668027. [8] Agust´ın Salgado and Javier S´anchez. A temporal regularizer for large optical flow estimation. In IEEE International Conference on Image Processing ICIP, pages 1233–1236, 2006. http://dx.doi.org/10.1109/ICIP.2006.312548. [9] Javier S´anchez P´erez, Enric Meinhardt-Llopis, and Gabriele Facciolo. TV-L1 Optical Flow Estimation. Image Processing On Line, 2013:137–150, 2013. http://dx.doi.org/10.5201/ ipol.2013.26. [10] Andreas Wedel, Daniel Cremers, Thomas Pock, and Horst Bischof. Structureand motionadaptive regularization for high accuracy optic flow. In IEEE International Conference on Computer Vision, pages 1663–1668, September 2009. [11] C. Zach, T. Pock, and H. Bischof. A Duality Based Approach for Realtime TV-L1 Optical Flow. In Fred A. Hamprecht, Christoph Schn¨orr, and Bernd J¨ahne, editors, Pattern Recognition, volume 4713 of Lecture Notes in Computer Science, chapter 22, pages 214–223. Springer Berlin Heidelberg, Berlin, Heidelberg, 2007. 260