IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

1

Motion Detail Preserving Optical Flow
Estimation
Li Xu, Member, IEEE, Jiaya Jia, Senior Member, IEEE, Yasuyuki Matsushita, Senior Member, IEEE
Abstract—A common problem of optical flow estimation in the multi-scale variational framework is that fine motion structures
cannot always be correctly estimated, especially for regions with significant and abrupt displacement variation. A novel extended
coarse-to-fine (EC2F) refinement framework is introduced in this paper to address this issue, which reduces the reliance of flow
estimates on their initial values propagated from the coarse level and enables recovering many motion details in each scale.
The contribution of this paper also includes adaptation of the objective function to handle outliers and development of a new
optimization procedure. The effectiveness of our algorithm is demonstrated using the Middlebury optical flow benchmark and by
experiments on challenging examples that involve large-displacement motion.

!

1

I NTRODUCTION

The variational framework [18], together with coarse-tofine refinement [2], [23], is widely used in optical flow
estimation [10], [12]. In the Middlebury optical flow evaluation website [3], [4], almost all top-ranked methods adopt
this scheme.
However, the conventional coarse-to-fine warping framework has a fundamental limitation in handling motion
details. Brox et al. [9], in computing large-displacement
optical flow, pointed out that if flow structure is smaller
than its displacement, the latter may not be well estimated.
In this paper, we show that this issue can be even more
serious, as it also applies to small-displacement motion.
Taking Fig. 1 as an example, due to the camera motion, theforeground toy deer has motion significantly different fromthat of the background (average displacements d = −2 andd = 21 respectively). This example is very challenging forthe coarse-to-fine variational optical flow estimation.
As shown in Fig. 1(e), in a coarse level, the narrow neckdoes not exist and only the significant background motion isestimated. This makes the actual motion of the foregroundpixels in the finer scale (Fig. 1(f)) drastically different
from the initial estimate from the background, violating thelinearization assumption and accordingly leading to highly
unstable motion estimation. The final flow result shown in(c) includes considerable errors. This example discloses one
problem of the general coarse-to-fine variational model –that is, the inclination to diminish small motion structures
when spatially significant and abrupt change of the displacement exists.
We address this problem in this paper and propose
a unified framework for high-quality flow estimation in• L. Xu and J. Jia are with the Department of Computer Science andEngineering, The Chinese University of Hong Kong, Hong Kong.
E-mail: xuli@cse.cuhk.edu.hk, leojia@cse.cuhk.edu.hk• Y. Matsushita is with Microsoft Research Asia, Beijing, China.
Email: yasumat@microsoft.com

Coarse Level

(a)

(b)
(e)

Finer Level
(c)

(d)

(f)
(g)

Fig. 1. Motion detail preserving problem. (a)-(b) Two input
patches. (c) Flow estimate using the coarse-to-fine variational setting. (d) Our flow estimate. (e)-(f) Two consecutive
levels in the pyramid. Flow fields are visualized using thecolor code in (g). The input data is from [29].

both large and small displacement settings. Central to our
method is a novel selection scheme to compute extended
initial flow vectors in each image level. This makes thefollowing optimization not completely rely on the resultat the previous scale, and is thus capable of refining
the estimation correctly in a top-down fashion. Our flow
result shown in Fig. 1(d) contains small structures. More
examples are included in Section 5.
This paper also contributes in the following ways. First,
we use robust sparse features, together with patch matching,
to produce extended flow initialization, which helps enforce
the linearization condition in the variational setting. Second,in the flow estimation model, we propose the selective
combination of color and gradient constraints in defining
the data term, robust to outliers. Third, we propose a fastvariable-splitting-based optimization method to refine flow
maps. It is highly parallel.
Finally, we employ the Mean Field approximation to
enable solving the objective function, which involves bothIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

discrete and continuous variables, commonly regarded as
challenging to handle. Extensive experiments visually andquantitatively validate the performance of our approach forboth large- and small-displacement motion.
This manuscript extends its conference version [42] withthe following major differences. (1) We provide morediscussion and derivation details of the data term and its
Mean-Field approximation. (2) The extended coarse-to-finescheme is further generalized in order to handle non-rigid
motion with dense patch matching. (3) The occlusion isprogressively handled in each scale. (4) We have experimented with more challenging examples in this paper.
The rest of the paper is organized as follows. Section 2
reviews related work. In Section 3, we introduce the flowenergy. An extended coarse-to-fine framework together withthe efficient solver is presented in Section 4. Section 5shows results for both large- and small-displacement opticalflow. We conclude this paper in Section 6.

2

R ELATED WORK

Following the variational model of Horn and Schunck [18],modern optical flow estimation is usually posed as anenergy minimization problem, with the energy function
containing a data term and a smoothness term.
One important improvement over the original variationalmodel is the introduction of robust statistics for both theenergy terms. Black and Anandan [7] replaced the quadratic
penalty functions in [18] by non-convex robust ones to reject outliers. 1 -norm, or its variation (e.g., the Charbonnierfunction), is also commonly used [10], [12], [37], [43].
Learning-based methods construct the distribution fromempirical data. In particular, Roth and Black [26] learned
spatial smoothness using Field-of-Experts, and combined itwith a Charbonnier data term; Sun et al. [34] proposeda learning framework to fit distributions using Gaussian
Scale Mixture (GSM). In [33], Sun et al. empiricallydemonstrated that the simple Charbonnier function (1 norm) actually outperforms other highly non-convex robust
functions due to its convex property.
Efforts also have been put into improving the opticalflow constraints. Haussecker and Fleet [17] proposed aphysical constraint to model brightness change. Wedel etal. [37] proposed a structure-texture decomposition method
to reduce the discrepancy between two frames caused byillumination change. Lempitsky et al. [20] computed thematching cost only using high frequency components. Prefiltering on the input images was suggested in [34] and [25]to handle illumination variation. These models are flexible,
but at the same time require pre-processing or advanced
optimization to solve complex objective functions.
In [10], Brox et al. introduced the gradient constancy
constraint to complement the brightness constraint. In [12],separate penalties are imposed on the brightness and gradient constraints. Zimmer et al. [45] further employed thenormalized brightness and gradient constraints. We willshow later that the way to combine the brightness andgradient terms can be improved by a selection model.

2

To preserve motion discontinuities, anisotropic- [39],steerable- [34], [44], [45], and adaptive-smoothness [1],[36] terms were studied. Segmentation information wasincorporated to regularize flow estimates in [19], [24],[41]. Recently, the non-local smoothing strategy [33], [38]demonstrated the potential to handle displacement discontinuities and occlusion, which is tightly linked to explicit
refinement of the flow field using image filtering [28], [40].
Almost all the above methods rely on coarse-to-fine
warping to deal with motion larger than one pixel [2], [6].
As discussed in Section 1, this strategy could fail to recoversmall-scale structures. Handling incorrect initialization by
adapting windows for stereo matching [31] is a solution.
It however assumes at least that the nearby disparities arecorrectly initialized, which may not be true for small-scale
structures that are totally eliminated in the coarse level.
Using discrete optimization, Lempitsky et al. [20] proposed fusing flow proposals obtained from different flow
estimation methods with various parameter settings. It iseffective at finding the optimal values among the given
proposals. But the sufficiency and optimality of the proposals cannot be controlled. Because the proposals are still
generated by the conventional coarse-to-fine warping, it ispossible that none of the proposals preserve small-scale
motion structure. In comparison, our method computes highconfidence flow candidates in each level, and thus is notentirely dependent on flow obtained from the coarse scale.
Related work also includes recent large-displacement
optical flow estimation [9], [11], where region-based descriptor matching was introduced. It is an effective method
except for occasional vulnerability to matching outliers, due
to the data term. As discussed in [11], descriptor matching
could decrease the performance in small-motion regions.
By extending the numerical scheme of [43] and by
searching possible values to minimize the data energy [32],large displacement optical flow estimation can be achieved.
As the smoothness prior is not enforced, the results can possibly be noisy and lack sub-pixel accuracy. In this paper, an
extended coarse-to-fine method is proposed, which can significantly improve both the large- and small-displacement
optical flow estimation in a unified framework.

3

O PTICAL F LOW M ODEL

We introduce in this section our objective function. We baseour data penalty function on the 1 norm to reject outliers
and use the Total Variation (TV) for regularization.
3.1 Robust Data Function
As the color constancy constraint is often violated whenillumination or exposure changes, combining the gradientconstraint was adopted [10], [12]. Denoting by u = (u, v)T
the flow field that represents the displacement between
frames I1 and I2 , one choice of the data term for flow
estimation is
1
ED (u) =I2 (x + u) − I1 (x) +2
x1
τ ∇I2 (x + u) − ∇I1 (x),2(1)

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE
10

α=1
α=0
α=0.5
eff. energy

8
6

10
6

4

4

2

2

discontinuity [1], [36]:

α=1
α=0
α=0.5
eff. energy8ω(x) = exp(−∇I1 κ ),0
0-5 -4 -3 -2 -1 0 1 2 3 4 5 -5 -4 -3 -2 -1 0 1 2 3 4 5(a)(b) Data costs for P1(c) Data costs for P2

in the “RubberWhale” example, where two points P1 and P2
are highlighted. (b)-(c) Plots of different data costs (heightsof the points) for P1 and P2. The ground truth displacement
is moved to 0 in the horizontal axis for ease of illustration.
2

where x ∈ Z indexes the 2D coordinates, τ is a weightbalancing the two matching costs. ∇ is the discrete approximation of the gradient operator. This function, due tothe addition of two terms, is less accurate in modeling pixel
correspondence than only using one out of the two terms.
Fig. 2 shows an example where the patch in (a) contains
two points P1 and P2. Their data cost distributions withrespect to different displacement values are plotted in (b)and (c) respectively (ground truth displacements are shiftedto 0). It is noticeable that the color constraint (blue curvein (b)) does not produce the minimum energy near the
ground truth value because the color constancy is violated
given point P1 moving out of the shadow. Adding thecolor and gradient terms using Eq. (1) also results in an
undesirable distribution (dashed magenta curve) as the costat the ground truth point is not even a local minimum.
Similarly, in Fig. 2(c), only the color constancy holds
because point P2 undergoes rotational motion, which alters
image gradients. It is not ideal as well to add the twoconstraints in the data function definition.
The above analysis indicates that a good model should
only use the more fitting constraint, but not both ofthem. We accordingly define a binary weight map
α(x) : Z2 → {0, 1} to switch between the two terms. The
new data function is expressed as
α(x)I2 (x + u) − I1 (x) +ED (u, α) =x

(2)When α(x) = 1, the gradient constraint is favored. Otherwise, we select color constancy. Our empirical investigation
provided in Section 5 shows that this model can produce
higher quality results than various alternatives.
3.2 Edge-Preserving Regularization
The regularization term for optical flow estimation is generally designed to be edge preserving [34], [36], [45]. We
define our smoothness term asω(x)∇u(x),(3)
ES (u) =x

where ∇u(x) is the common TV regularizer. ω(x) is
the simple structure adaptive map that maintains motion(4)

where we set κ = 0.8 in our experiments. For simplicity,
we use the brightness derivatives (2 channels) to compute∇I1 κ . The final objective function is defined asE(u, α) = ED (u, α) + λES (u),Fig. 2. Data cost distributions for two points. (a) A patch(1 − α(x))τ ∇I2 (x + u) − ∇I1 (x).

3(5)where λ is the regularization weight.
3.3

Mean Field Approximation

Minimizing Eq. (5) involves simultaneously computing two
fields: continuous u and binary α, which is computationallychallenging. We employ the Mean Field (MF) approximation [14] to simplify the problem by first canceling out thebinary process by integration over α. The probability of aparticular state of the system is given by
1 −βE(u,α)
e
,
(6)
Z
where β is the inverse temperature and Z is the partition
function, defined as e−βE(u,α) .
(7)Z=P (u, α) ={u} {α=0,1}

We then compute the sum over all possible αs with thesaddle point approximation (see Appendix in the supplementary file for the derivation), yielding:1ln(e−βDI (u,x) + e−βD∇I(u,x) ), (8)E ef f (u) = λES (u) −
β
x
and
where
DI (u, x) = I2 (x + u) − I1 (x)D∇I (u, x) = τ ∇I2 (x + u) − ∇I1 (x).
It
indicatesthat the flow estimate by minimizing Eq. (8) is actually the
Mean Field (MF) approximation of minimizing Eq. (5).
The effective energy is therefore written asef f(u) + λES (u),E ef f (u) = ED

(9)where the effective data function is
 1
ef f
ED
(u) =
− ln(e−βDI (u,x) + e−βD∇I(u,x) ). (10)βxThe optimality of Eq. (9) does not depend on the estimateof α. Eq. (10) defines a robust function and β plays a keyrole in shaping it. When β → 0, Eq. (10) acts as the averageof the two data costs in Eq. (1), while β → ∞ leads to thelower envelope of the two costs in Eq. (10). We show in
Fig. 3 several examples on how the effective function isaffected by varying β. Fig. 3(a) contains plots with different
α values, same as the ones shown in (b). In (b)-(d), weshow distributions of the effective data costs by varying β.
Note that a small β makes the distribution (plotted in (b))
close to the original one with α = 0.5 (shown in (a)) whilea relatively large β (shown in (d)) yields the distribution
approaching the lower envelope of the costs with α = 0 and
α = 1, which is what we need for accurate flow estimation.
The effective data costs (with β = 5) are also plotted inIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

α=1
α=0
α=0.5

α(x) (see Appendix in the supplementary file), and thus
can be updated once u is obtained. It has an effect similarto α(x) (given in Eq. (2)) in constraint selection.

β=0.01

4
0

5

0

5

(b) β = 0.01(a) Original plots

β=5

β=0.3

0

5

0

(c) β = 0.345(d) β = 5

Fig. 3. Effective data cost distributions with different β
values. (a) shows data costs with different α values. (b)-(d)show effective data costs by varying β.
Input: a pair of images for optical flow estimation
1. Construct pyramids for both of the images and set the initiallevel l = 0 and ul = 0 for all pixels.
2. Propagate ul to level l + 1.
3. Extended Flow Initialization (Section 4.1)
3.1. Detect and match SIFT features in level l + 1.
3.2. Perform patch matching in level l + 1.
3.3. Generate multiple flow vectors as candidates.
3.4. Optimize flow using QPBO (Eq. (9)).
4. Continuous Flow Optimization (Section 4.2)
4.1. Compute the ᾱ map (Eq. (12)).
4.2. Solve the TV/1 energy function in Eq. (15).
5. Occlusion-aware Refinement (Section 4.3)
6. If l = n − 1 where n is the total number of levels, l = l + 1
and go to Step 2.
Output: the optical flow field

TABLE 1O PTIMIZATION F RAMEWORK

Traditional optical flow estimation, ascribed to the use of
the variational setting, relies excessively on the coarse-tofine refinement. As discussed in Section 1, this process
could fail to recover ubiquitous fine motion details given
the possible large discrepancy between the initial flow andthe ground truth displacements in each image level.
In this section, based on E ef f and ᾱ, we propose an
iterative method to optimize Eq. (5). Specifically, because
ef f(u) is independent of ᾱ, we first infer multiple highED
confidence flow candidates and apply discrete optimization
to select the optimal ones. With this result, ᾱ in Eq. (12) canthen be quickly estimated. We finally improve the subpixel
accuracy of flow with the estimated ᾱ using continuous
optimization. This procedure is found to be surprisingly
effective in dampening estimation errors.
Our overall algorithm is sketched in Table 1 based on
iteratively processing images in a top-down fashion. The
steps are detailed further below.
4.1

Extended Flow Initialization

We address the general flow initialization problem in each
image scale by finding multiple extended displacements(denoted as {uv0 , ..., uvn }) through sparse feature matching and dense patch matching to improve estimation in
uc , which is the flow field computed in the immediately
coarser level. The following steps are adopted to obtain theextended displacements.

Method Overview

Fig. 2(b) and (c) using the green crossed curves. They
are coincident with the smallest-value curves. Our method
always keeps β ≥ 1 empirically.
We optimize Eq. (9) using an iteratively reweighted
optimization strategy. The difficulty of minimizing Eq. (9)stems from the non-convex data function. Taking the partial
derivative with respect to the variable u yields

ef f
(u) =
ᾱ(x)∂u DI + (1 − ᾱ(x))∂u D∇I , (11)∂u EDxwhere ᾱ(x) is the flow-dependent weight, written asᾱ(x) =

1.
(12)
1 + eβ(DI (u,x)−D∇I (u,x))
It indicates that the energy can be minimized by iteratively
updating ᾱ in the outer loop and by solving for u with thecomputed weights afterwards. In addition, although Eq. (9)is non-convex and difficult to solve using continuous optimization, there is no obstacle to apply discrete optimization
if candidate labels can be obtained. We propose a robustalgorithm, described in the next section, to estimate u.
The solver can also be interpreted from another perspective. Note that ᾱ(x) is actually the MF-approximation of

4.1.1 SIFT feature detection
SIFT feature detection and matching [22] can efficiently
capture large motion for objects undergoing translationaland rotational motion. Instead of computing a dense descriptor field as in scene matching [21], we only employ
sparse matching of discriminative points, which avoids
introducing many ambiguous correspondences and outliers.
One example is shown in Fig. 4(b). Note that some matches
could still be wrong. But this is not a problem as we will
eventually employ discrete optimization to only select themost credible candidates.
4.1.2 Selection
The computed displacement vectors by feature matching
are denoted as {s0 , ..., sn }, as shown in Fig. 4(b). They
are new potential flow candidates except those that alreadyexist in the initial flow field uc (shown in Fig. 4(c)). To
robustly screen out the duplicated vectors, we computethe Euclidean distance between each si and all ucj s where
pixel j is within a 5 × 5 window centered at the reference
feature si . If all results are greater than 1 (pixel), we regard
si as a new flow candidate. We repeat this process for
all is, and denote the m remaining candidate vectors as{sk0 , ..., skm−1 }, as shown in Fig. 4(d).

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

5

(g)

(c)

(a) u0 map

(b)

(e)
(d)

(b) ᾱ(x) map(c) ur map(h)(d) Close-up (e) Close-up
(a)

(f)

(i)

Fig. 4. Extended flow initialization. (a) Two input frames. (b)One of the images overlaid with the computed feature motion
vectors si . (c) Flow field uc propagated from the coarse level.
(d) New displacements {sk0 , ..., skm−1 } computed using (b)and (c). (e) New displacement maps. Each usi is expandedfrom ski and therefore is a constant-value map. (f) Dense
nearest-neighbor patch matching field un . (g) Optimized flow
map u0 with respect to all candidates in the current image
scale. (h)-(i) show close-ups of (c) and (g).

This strategy significantly reduces the system dependence on the coarse-scale flow estimation. It is notable aswell that feature matching initially produces many vectors
distributed in the whole image, as shown in Fig. 4(a);but they reduce to less than 15 candidates after local
comparison with uc in the given example. Only the mostdistinct flow vectors are retained.
4.1.3

Expansion

The m remaining vectors {sk0 , ..., skm−1 } represent possible missing motion in the present flow field uc . To
determine whether or not they are better estimates to replace
the original ones in uc , we expand each displacement vector
ski to a constant-value flow field usi for further fusion. The
fields are shown in Fig. 4(e).
4.1.4Patch matching

SIFT Feature matching, albeit very effective, sometimes
still misses motion vectors. It is because small texturelessobjects may not have distinct features, making their detection problematic. Another main reason is that to let SIFT
descriptors gather enough information for 128-dimension
feature vector formation, the patches on which they operateshould at least contain 16 × 16 samples as suggested.
The size could be too large for non-rigid motion as edge
statistics may change a lot for successive two frames.
We resort to dense nearest-neighbor patch matching foramelioration. The patches we use can be as small as
5 × 5. They are more flexible to describe motion of small
textureless regions, as shown in Fig. 4(f). Specifically, we
compute the matching field un by minimizing energy E(un , x) =
Ik2 (y + un (x)) − Ik1 (y)2 , (13)y∈N (x) k

Fig. 5. Continuous optimization. Errors are further reducedin this step. (d) and (e) show close-ups of (a) and (c).

where Ik ∈ {Ir , Ig , Ib , ∂x I, ∂y I}, denoting a total of 5 colorand gradient channels. N (x) is a 5 × 5 window centeredat x. Although noise is generated by this method, it can bequickly rejected in the following optimization step with thecollection of a set of flow candidates for each pixel.
The energy (13) was employed in [13] as well. But
linearization was performed eventually in [13], confining
only local refinement. In comparison, we do not impose
any smoothness constraint at this stage. So estimates forvery large displacement can be obtained.
4.1.5 Matching field fusion
The m + 1 new motion fields {us0 , ..., usm−1 , un }, together
with the original uc , comprise several motion candidates
for each pixel in the present image scale. Selection of the
optimal flow among the m+2 candidates for each pixel is alabeling problem, with the objective function in Eq. (9). It
can be solved by discrete optimization efficiently because
on the one hand the number of candidates is small, thanks
to the carefully designed selection process; on the other
hand, Eq. (9) does not involve α, simplifying computation.
We adopt the Quadratic Pseudo-Boolean Optimization(QPBO) [27] to solve this problem. The fusion movestep [20] is used to repeatedly fuse the candidates until
each gets visited twice. Also, to suppress the checker-boardlike artifacts commonly produced near motion boundaries
in discrete optimization, we employ the anisotropic representation of the TV regularizer ∇u = ∇u1 + ∇v1
with 8-neighbor discretization [15]. This method turns thechecker-board-like boundaries to octagons, a better approximation of the original smooth boundaries. The output is theflow map denoted as u0 . One result is shown in Fig. 4(g),which contains better recovered motion structure comparedto the field uc in (c). Close-ups are shown in (h) and (i).
Note that an alternative is to directly discretize theoriginal 2D solution space and fuse all candidate flows. It
however may suffer from expensive and possibly unstablecomputation because hundreds of labels can be produced
simultaneously in the original resolution.
4.2

Continuous Flow Optimization

We now refine flow u0 through continuous optimization,by iteratively updating ᾱ in Eq. (12) and u. The initial

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

flow field is taken into Eq. (12) to estimate ᾱ, as shown in
Fig. 5(b). Considering that Eq. (10) is highly non-convex,
we then take ᾱ back to Eq. (5) for optimization in thevariational model.
As color images are used, we still denote byIk ∈ {Ir , Ig , Ib , ∂x I, ∂y I} the set of channels included in thedata term and use αk ∈ {ᾱ, ᾱ, ᾱ, (1 − ᾱ)τ, (1 − ᾱ)τ } torepresent the corresponding weights. The energy in Eq. (5)is thus written asE(u) =αk (x)Ik2 (x + u) −
(14)
x

k
k
I1 (x) + λ(x)∇u(x),

x

k

(15)

given small du. In Eq. (15),Ix
Iy= ∂x I2 (x + u0 ),= ∂y I2 (x + u0 ),It= I2 (x + u0 ) − I1 (x).

To preserve motion discontinuity, we employ the rotational
invariant isotropic form of the TV regularizer, written as∇u = (∂x u)2 + (∂y u)2 + (∂x v)2 + (∂y v)2 . (16)
Our Solver We propose decomposing the optimization
into three simpler problems, each of which can have theglobally optimal solution. The key technique is a variablesplitting method [35] with auxiliary variables p and w,
representing the substituted data cost and flow derivatives respectively, to move a few terms out of the nondifferentiable 1 -norm expression. This scheme is found tobe efficient and is crucial to produce high quality results.
The derivatives of each flow vector comprise four elements, i.e.,∇du = (∂x du, ∂y du, ∂x dv, ∂y dv)T .
For each element, we introduce a corresponding auxiliary
variable. The set of the variables is denoted as
w = (wdux , wduy , wdvx , wdvy )T .
Then Eq. (15) is transformed to 1Ikx du + Iky dv + Ikt − pk 2 + αk pk  +2ηxkInput: images Ik , initial flow field u0 , weights αk .
Perform linearization at u0η ← η0repeatCompute pk using Eq. (19)θ ← θ0repeatCompute w using Eq. (21).
Compute du by solving Eq. (22).
θ ← θ/3
until θ < θmin
η ← η/3
until η < ηminur = u0 + du
Output: refined flow field ur .

TABLE 2where λ(x) := λω(x). With the initial flow u0 estimated in the previous step, we solve for the increments
du = (du, dv)T by minimizing Eq. (15). The final flow
vector is u = u0 + du. By convention, the Taylor expansion of Eq. (15) at point x + u0 yields
E(u) =αk (x)Ikx du + Iky dv + Ikt  +λ(x)∇(u0 + du)(x),61∇du − w2 + λ∇u0 + w. (17)2θ1
In this function, 2ηIkx du + Iky dv + Ikt − pk 2 + αk pk 
kencourages p
to approach Ikx du + Iky dv + Ikt , and0122θ ∇du − w + λ∇u + w2 makes w similar to

Algorithm for continuous flow optimization∇du. It can be observed as well that Eq. (17) approaches
Eq. (15) when θ → 0 and η → 0. Our algorithm proceedswith the following iterations with initial u := u0 .
1. Fix u to estimate p. The simplified objective function is 1minIkx du + Iky dv + Ikt − pk 2 + αk pk . (18)2ηxkSingle variable optimization can be used in this step. The
optimal solution is given by the shrinkage formula [16]pk = sign(ok ) max(|ok | − ηαk , 0),(19)where ok := Ikx du + Iky dv + Ikt is the flow constraint.
2. Fix u to estimate w. The function reduces to 1min∇du − w2 + λ(x)∇u0 + w2 .
2θ
x

(20)

Similarly, the following solution can be obtained by theshrinkage formula
wdux = max(∇u2 − θλ, 0)∂x u− ∂x u0 ,
∇u2

(21)

where u = u0 + du. Solutions for wduy , wdvx , and wdvy
can similarly be derived. The computation in this step isalso quick and is highly parallel by nature.
3. Fix w, p and solve for u. The objective function is 11Ikx du+Iky dv +Ikt −pk 2 + ∇du−w2 .
min
2η
2θ
x
k
(22)
It is quadratic and the corresponding Euler-Lagrange equations of Eq. (22) are linear w.r.t. du and dv. A globally
optimal solution can be obtained by solving the linear
system in this step.
Our method iterates among optimizing (19), (21), and(22) until convergence. Note that cost function decomposition with auxiliary variables was used in [38], [39], [43] forflow estimation. Their steps use the primal-dual solvers. In
comparison, our scheme consists of a set of simpler subproblems, each with guaranteed global optimality. It thus

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

7(a) “Hydrangea”(b) n=1, AAE=24.37 (c) n=2, AAE=2.52 (d) n=5, AAE=2.45(e) ground truth(f) n=2, AAE=4.83 (g) n=5, AAE=2.63 (h) n=10, AAE=2.60(i) Energy plot

Fig. 6. Continuation scheme. (b)-(d) show our results obtained using the algorithm in Table 2 with the continuation scheme.
n is set to 1, 2, and, 5 respectively. The error is already very small when n = 2. (f)-(h) show results with fixed η = 0.1 andθ = 0.01 in all iterations. AAE stands for “average angular error” [3]. (i) Energy decreasing w.r.t. the number of iterations withand without continuation.
differs from previous methods in the way of formulating the
problem and of proposing the solver to each sub-problem.
In practice, θ and η are critical parameters that should besmall. It was found that fixing them to constants typicallyresults in slow convergence. We thus adopt the continuation
scheme [16] for speedup, which initially sets θ and η tolarge values to allow warm-starting, and then decreases
them in iterations toward the desired convergence. Ouralgorithm is sketched in Table 2, where ηmin and θminare set to 0.1 and 0.01 respectively. η0 and θ0 are therespective initial values, configured as η0 = 3n × ηmin andθ0 = 3n × θmin , where n controls the number of iterations.
Fig. 5(d) and (e) show flow fields before and after thecontinuous refinement in an image scale. We denote byur the refined flow field.
Fig. 6 demonstrates the effectiveness of this continuation
scheme (that is, by altering η and θ in iterations) andcompares results obtained with and without using it. We
set different iteration numbers in experiments. The top row
shows results with n = 1, n = 2, and n = 5 using thecontinuation scheme. The bottom row contains estimatesusing the algorithm shown in Table 2, by fixing η = 0.1
and θ = 0.01 in all iterations. Energy decreasing w.r.t. thenumber of iteration is plotted in (i). It is clear from thecomparison that our algorithm with the continuation scheme
converges more efficiently.

(a) Occlusion and flow(b) urFig. 7. Occlusion-aware refinement. (a) Flow estimate
overlaid with the occlusion map (o(x) > 0.5). (b) and (c) show
results before and after the final refinement in an image
scale.
where f(x + u(x)) is the count of reference pixels mapped
to position x + u(x) in the target view using forward
warping. Tl,h (a) is a function that truncates the value ofa if it is out of the range [l, h]. Eq. (23) indicates if there
exist more than one reference pixel mapping to x + u(x),the occlusion label for the reference x is set. Although thissimple method sometimes fattens the occlusion region, it
seldom leaves out true occluded pixels, and thus is usefulin the final flow estimation. In practice, we apply a small
Gaussian filter on the computed o(x) to reduce noise.
Our measure of the data confidence based on the occlusion detection is expressed asc(x) = max(1 − o(x), 0.01).

4.3 Occlusion-Aware Refinement
Motion vectors for occluded pixels generally cannot be
determined due to the lack of correspondences. In this step,
we handle occlusion in the computed flow field. Althoughcross-checking is effective in occlusion detection, it needs
to compute optical flow bidirectionally. Our strategy isbased on an observation that multiple pixels mapping to
the same point in the target image using forward warpingare possibly occluded by each other.
Thus, we detect occlusion using the mapping uniquenesscriterion [8], expressed aso(x) = T0,1 (f(x + u(x)) − 1),(23)(c) Refined(24)

The value 0.01 is to make c(x) always larger than 0.
The metric is used in the following two ways to improveflow estimation in the occluded regions. First, we explicitly
perform cross bilateral filtering for the detected occluded
pixels where o(x) > 0.5. Each pixel is further weighted
by the measure c(x) so that occluded pixels have weakerinfluence in filtering. This scheme was shown to be effectivein occlusion handling [28], [40] and was used in defining
the flow function [33], [38].
Second, based on the fact that we should not trust thedata term with large o(x), the energy function is updatedwith respect to the occlusion confidence, which makes flow
computation for the occluded pixels depend more on theIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

(a)

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(b)

Fig. 8. Flow estimation error comparison. (a) Estimation
errors w.r.t. α. α = 1 and α = 0 indicate respectively
that only the color or gradient constancy constraint is used.
α = 0.5 refers to weighted addition of the two constraints. (b)Estimation errors w.r.t. β on the “RubberWhale” example.
local smoothness constraint:E  (u) = c(x)ED (u) + λES (u).

(25)

It can be efficiently optimized also with our solver.
The occlusion-aware flow refinement is applied at eachscale with the computed vectors from the continuous estimation step. The final result of the “Grove” example in
one image scale is shown in Fig. 7 where the detected
occlusion map is overlaid on the flow estimate. We compare
the ur maps obtained before and after our occlusion-aware
refinement in (b) and (c).

58Fig. 9. Visual comparison with different α settings. (a)and (b) show two image patches. (c) and (d) show flow
results computed using the color and gradient constraints
respectively. (e) is the ground truth flow field. (f) shows theresult with α = 0.5. (g) is the flow map obtained using ourselective combination model. (h) shows the ᾱ map.

(a) Frame 2(b) No OH(c) With OH(d) Frame 1(e) No OH(f) With OH

E VALUATION AND E XPERIMENTS

In this section, we present our results in both small- andlarge-displacement settings. τ in Eq. (2) is set to 1/1.4
to normalize the color and gradient constraints, which islearned from the Middlebury training image set by setting
the color and gradient costs to be equal. In order to reduce
the sampling artifacts in Eq. (12), we filter DI and D∇Iwith a small Gaussian kernel with the standard deviation
1.0. β, λ, η, and θ are empirically set to 5, 12, 0.1,
and 0.01 respectively. For feature detection, we use theimplementation of Lowe [22]. Matches are retained only
if ratios between the best and the second best matchingscores are smaller than 0.6. For patch matching, we adopt
the randomized nearest-neighbor method approximation [5]with patch size 5 × 5.
5.1 Evaluation of the Data Term
We evaluate the selective combination strategy in defining
the data cost function. We compare our method with thoseusing fixed weights α = 0.5, α = 1, and α = 0 on
the Middlebury training set [4], where the ground truth
data are available. To demonstrate the alpha influence not
involving our other steps, we employ the classic coarse-tofine warping framework. The errors are listed in Fig. 8(a),calculated on the two representative examples “RubberWhale” and “Urban2”. It can be noticed that the averageangular error (AAE) for “Urban2” is small when using thecolor constraint alone while the gradient constraint is morefavored in “RubberWhale” due primarily to illumination
variation. Simply adding these two constraints (α = 0.5)

Fig. 10. Flow estimation of a car moving away from thecamera. (a) and (d) are two input images. (b) and (e) showrespectively the flow field produced without occlusion handling and the backward warping result. (c) and (f) are theflow field and the warping result with our occlusion handling.
“OH” stands for occlusion handling.

produces AAE in between. Our method locally selects themore optimal term and thus performs better.
Fig. 8(b) shows how the estimate changes with respect
to different β for the “RubberWhale” example. β = 0
corresponds to weighted addition of the two normalized
constraints. Note that the average error of the flow field
decreases quickly with the increase of β, in line with ourunderstanding (explained in Section 3.3).
In Fig. 9, we show a visual comparison. Red arrows in(a) and (f) indicate pixels violating the color constancy
assumption. The blue arrows highlight the edge of thewheel, of which the gradient varies. (c) and (d) show resultsby respectively setting α = 1 and α = 0. (f) shows theresult with α = 0.5, where problems caused by using eitherof the constraints is still present. Our selective combination
model helps robustly reject outliers, as shown in (g).
For quantitative comparison, a series of experiments withdifferent optimization strategies are conducted, varying
from traditional coarse-to-fine to our full optimization with
EC2F and occlusion refinement. The error statistics arelisted in Table 3, where “F” represents setting α = 0.5 andIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

Method
F+C2F
A+C2F
F+EC2F
A+EC2F
F+EC2F+o
A+EC2F+o
F+EC2F+O
A+EC2F+O

RubberWhale
3.29
2.95
2.95
2.61
2.93
2.59
2.92
2.59

Hydrangea
2.25
2.02
2.13
2.06
2.08
2.02
2.08
1.97

Dimetrodon
2.58
2.50
2.54
2.51
2.55
2.52
2.54
2.51

9

Grove2
2.26
2.09
2.07
2.00
2.04
1.92
1.97
1.88

Grove3
5.88
5.51
5.22
5.17
5.13
4.87
4.88
4.77

Urban2
3.61
2.81
2.66
2.40
2.46
2.22
2.40
2.15

Urban3
5.23
4.28
4.09
3.32
4.02
3.21
4.03
3.15

Venus
6.45
6.07
4.27
3.86
4.07
3.56
3.70
3.55

TABLE 3AAEs yielded by different strategies on the Middlebury optical flow training data

Initial Flow
C2F
C+SIFT
C+PM
All

RubberWhale
2.71
2.61
2.66
2.59

Hydrangea
2.00
1.97
2.00
1.97

Dimetrodon
2.51
2.50
2.51
2.51

Grove2
2.08
1.94
1.95
1.88

Grove3
5.19
4.79
4.83
4.77

Urban2
2.77
2.17
2.27
2.15

Urban3
4.14
3.20
3.87
3.15

Venus
5.87
3.65
4.13
3.55

TABLE 4AAEs on the Middlebury training data under different flow initialization

Fig. 11. The average end-point errors (EPEs) on the benchmark data as of Oct, 2010, copied from the Middlebury
website [4]. Our method is denoted as “MDP-Flow2”.

“A” stands for our adaptive α scheme. “C2F” and “EC2F”
represent the classic and extended coarse-to-fine schemes
respectively. As described in Section 4.1, “EC2F” uses
extended flow initialization at each scale. Our method yields
consistent quality improvement over other alternatives.
5.2 Evaluation of Occlusion Handling
We also evaluate the occlusion-aware refinement step. The
bottom several rows of Table 3 list the statistics produced
without (“*+EC2F”) and with occlusion-aware refinement(“*+EC2F+O”). “*+EC2F+o” stands for occlusion handlingused in the early version of the system [42], where crossbilateral filtering is not employed. Both occlusion handling
methods yield reasonable results.
For the special case that an object moves away fromthe camera, the occlusion regions could be largely fattened
because multiple pixels on the object when it is near could
be mapped to one pixel when it is far. Even in this case,the flow estimates can still be refined because our occlusion
handling in essence seeks flow discontinuity alignment withimage edges.
We show in Fig. 10 an example. (b) and (e) show our flow
and backward warping results without occlusion handling.
The moving car is correctly reconstructed but the flow is

not accurate at the occluded region (rightmost part of thecar). With occlusion handling, the results are those shownin (c) and (f). The seemingly incorrect warping result in(f) in fact indicates correct handling of occlusion. The flownear the boundary is a bit noisy in (c), owing to the fattened
occlusion region jeopardizing proper flow regularization. It
is one of the limitations.
5.3

Evaluation of Extended Coarse-to-Fine

We also evaluate our coarse-to-fine framework with extended flow initialization. Specifically, we have tested 1)
classical flow initialization in the coarse-to-fine framework;2) flow initialization extended by SIFT feature matching
only; 3) extended flow initialization with patch matching; 4)
our flow initialization with both SIFT and patch matching.
Their abbreviations are “C2F”, “C+SIFT”, “C+PM” and“All”. The statistics are listed in Table 4. The results indicate that extended flow initialization (“C+SIFT”, “C+PM”,“All”) can greatly improve estimation.
5.4

Middlebury Optical Flow Benchmark

We now evaluate our method on the Middlebury opticalflow benchmark data. The table in Fig. 11 is copied in partIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

10(a) Input(b) Warping(c) Ground truth(d) Ours(e) LDOF [9](f) C2F [10](g) LDOF [9](h) [32](i) Ours(j) Ours

Fig. 12. Visual comparison on a large-displacement optical flow example from the HumanEva-II data set [30].

(a)

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Fig. 13. Visual comparison of the small-displacement optical flow results on two examples. (a) Our flow results. Closeups of (b) the input image, (c) ground truth flow, (d) our
estimate, and of results of (e) Brox et al. [10], (f) LDOF [11],(g) Zimmer et al. [45], (h) Werlberger et al. [38], and (i) of Sun
et al. [33], are shown.

from the evaluation website [4]. Our method, denoted as“MDP-Flow2”, ranked 1st at the time of submission (as ofOct. 29, 2010).
Regarding the running time, in our current CPU implementation, the whole program takes 420s to compute ahigh quality flow field for an image pair with resolution
640×480 in, for instance, the Urban sequence. The running
time is reported on a laptop computer containing an Intel
Core i7 CPU @2.13GHz and 2GB Memory.
We show our flow results for two examples in Fig. 13(a)-(b). Methods of Brox et al. [10] that uses TV/1 model
for flow estimation, the large-displacement optical flow
estimator [11] that incorporates descriptor matching in thedata term, and of three top-performing methods that provide motion-discontinuity-preserving regularization terms,
produce results shown in (e)-(i).
5.5 Large-Displacement Optical Flow Estimation
Our method by nature can deal with large-displacement
flow, without any modification of the framework. One
example from the HumanEva-II benchmark dataset [30]is shown in Fig. 12. It contains significant articulated
motion of a running person. The fast foot movement cannot
be estimated correctly in the conventional coarse-to-finescheme [10], as shown in (f). (b) shows the backward
warping result based on our dense flow estimate. The closeups are shown in (d). Our method successfully recovers theshape of the left foot. The pixels in the occluded region are
simply unknown for all optical flow estimation methods.
The flow magnitude maps are shown in the second row.
The maps in (g) and (h) are produced by two representativelarge-displacement optical flow methods.
The flexibility of our method is boosted by patch matching especially for non-rigid large-displacement motion estimation. Fig. 14(a) and (b) show two frames. The duck head
undergoes very large motion. So it is not surprising that
other optical flow methods based on the traditional coarseto-fine scheme [10], [33] cannot cope with it well. The
large-displacement methods using descriptor matching [11],[21], [42] produce results shown in Fig. 14(g)-(l). There
are also errors. The flow and warping results produced
by extended coarse-to-fine with only patch matching areIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

11(a) Frame 1(b) Frame 2(c) Brox et al. [10](d) Brox et al. [10](e) Sun et al. [33](f) Sun et al. [33](g) SIFT Flow [21](h) SIFT Flow [21](i) LDOF [11](j) LDOF [11](k) MDP-Flow [42](l) MDP-Flow [42](m) Patch Matching(n) Patch Matching(o) Ours(p) Ours

Fig. 14. A challenging example for large-displacement optical flow estimation. (a) and (b) show two input images. Seven
flow estimation and the corresponding backward warping results are shown in (c)-(p).

shown in (m) and (n). Although the field is noisy, it roughly
captures the head motion. Our final flow estimate, yieldedwith the complete EC2F scheme that involves both patch
matching and feature matching, is shown in (o). Its quality
is much higher.
Another example is shown in Figs. 15 and 16, which
is a low-frame-rate sequence containing a football player.
Fig. 15 contains the results of the conventional coarse-tofine warping method [10], the large-displacement estimator [11], and of our method. Fig. 16 shows a few resultsin the sequence. All examples demonstrate that in termsof handling large motion of small-size regions, our method
reduces the dependence on the linearization condition in
the variational model and thus can generate good results.

6

D ISCUSSION AND C ONCLUSION

We have presented a new optical flow estimation framework
to reduce the reliance on the coarse level estimation in thevariational setting for small-size salient motion estimation.
Differing from previous efforts mainly to improve themodel, we instead revise flow initialization in the coarse-tofine setting, which yields a unified framework to preserve
motion details in both small- and large- displacement
scenarios. The proposed method also takes advantage of theaccurate variational coarse-to-fine framework and of nonlocal search/matching. Other main contributions include theselective combination of the color and gradient constraints,
sparse feature matching and dense patch matching to collectIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

12(a) Frame 2(b) Brox et al. [10](c) LDOF [11](d) Ours(e) Frame 1(f) Brox et al. [10](g)LDOF [11](h) Ours

Fig. 15. Large-displacement optical flow results. (f)-(h) show the backward warping results based on the flow estimates in(b)-(d) respectively.

Fig. 16. Optical flow estimation in consecutive frames in a low-frame-rate sequence. First row: two-body-overlaid images tovisualize the large displacement. Second row: our flow estimates. Third row: magnitude maps.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

13(a) Frame 2(b) Ours(c) Ours+user input(d)(e) Frame 1(f) Ours(g) Ours+user input(h)Fig. 17. A challenging example. (a) and (e) are two input frames. (b) and (f) show our automatically computed flow field andthe backward warping result. (c) and (g) are the flow field and the backward warping result with simple user indication of two
additional pairs of correspondence, shown as the red and blue dots in (a) and (e). (d) and (h) are close-ups of (f) and (g).

appropriate motion candidates, the mean field approximation to simplify optimization, and a variable splitting
technique to enable fast and reliable flow estimation. Our
future work will be system acceleration using GPU.
Limitations There are several limitations. First, althoughsparse feature matching and dense patch matching complement each other in proposing new flow candidates, they
could still be insufficient especially for motion in textureless or regularly-patterned regions, where large matching
ambiguity could occur. Other information such as simple
user input may help.
We show one example in Fig. 17. (a) and (e) are input
frames where the boy’s right leg and arms undergo large
motion. (b) and (f) show our estimated flow field as wellas the backward warping result. Primary large-displacement
motion (of the shoe, for example) is correctly computed
except for the left arm. Note that textureless regions not
only fail feature detection, but also create ambiguities fornearest-neighbor matching. In this example, we manuallyspecify two corresponding points in the input frames (thered and blue dots in (a) and (e)), and then take thedisplacement vectors as new constant flow candidates toimprove estimation. Final results are shown in (c) and (g)with close-ups in (d) and (h). They indicate that simple
user interaction can decisively improve flow estimation inchallenging regions.
Second, motion inference for large occluded regions isstill an open problem due to lack of correspondence. Ourcurrent occlusion handling relies on a heuristic smoothness assumption, which could fail in texture- or color-richregions when occlusion is significant. Incorporating otherclues, such as color segmentation, may remedy the problem.

ACKNOWLEDGMENTS
The authors would like to thank the associate editor andall the anonymous reviewers for their time and effort. This
work is supported by a grant from the Research Grants
Council of the Hong Kong SAR (Project No. 412911) andby NSF of China (key project No. 61133009).

R EFERENCES[1]

L. Alvarez, J. Esclarin, M. Lefebure, and J. Sanchez. A pde model
for computing the optical flow. In Proc. XVI Congreso de Ecuaciones
Diferenciales y Aplicaciones, pages 1349–1356, 1999.
[2] P. Anandan. A computational framework and an algorithm for themeasurement of visual motion. International Journal on Computer
Vision (IJCV), 2:283–310, 1989.
[3] S. Baker, D. Scharstein, J. Lewis, S. Roth, M. J. Black, andR. Szeliski. A database and evaluation methodology for optical flow.
In ICCV, 2007.
[4] S. Baker, D. Scharstein, J. Lewis, S. Roth, M. J. Black, andR. Szeliski. A database and evaluation methodology for opticalflow. Technical report MSR-TR-2009-179, Microsoft Research,
2009. http://vision.middlebury.edu/flow/.
[5] C. Barnes, E. Shechtman, A. Finkelstein, and D. B. Goldman.
Patchmatch: a randomized correspondence algorithm for structural
image editing. ACM Transactions on Graphics (TOG), 28(3), 2009.
[6] J. R. Bergen, P. Anandan, K. J. Hanna, and R. Hingorani. Hierarchical model-based motion estimation. In ECCV, pages 237–252,
1992.
[7] M. J. Black and P. Anandan. The robust estimation of multiple
motions: Parametric and piecewise-smooth flow fields. Computer
Vision and Image Understanding (CVIU), 63(1):75–104, 1996.
[8] M. Z. Brown, D. Burschka, and G. D. Hager. Advances incomputational stereo. IEEE Transactions on Pattern Analysis andMachine Intelligence (PAMI), 25(8):993–1008, 2003.
[9] T. Brox, C. Bregler, and J. Malik. Large displacement optical flow.
In CVPR, 2009.
[10] T. Brox, A. Bruhn, N. Papenberg, and J. Weickert. High accuracy
optical flow estimation based on a theory for warping. In ECCV,
volume 4, pages 25–36, 2004.
[11] T. Brox and J. Malik. Large displacement optical flow: descriptor
matching in variational motion estimation. IEEE Transactions on
Pattern Analysis and Machine Intelligence (PAMI), 33(3):500–513,
2011.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE[12] A. Bruhn and J. Weickert. Towards ultimate motion estimation:
Combining highest accuracy with real-time performance. In ICCV,
pages 749–755, 2005.
[13] A. Bruhn, J. Weickert, and C. Schnörr. Lucas/kanade meetshorn/schunck: Combining local and global optic flow methods.
International Journal on Computer Vision (IJCV), 61(3):211–231,
2005.
[14] D. Geiger and F. Girosi. Parallel and deterministic algorithms formrfs surface reconstruction and integration. A.I. Memo 1114, MIT,1989.
[15] D. Goldfarb and W. Yin. Parametric maximum flow algorithms
for fast total variation minimization. Technical Report 07-09, Rice
University, 2007.
[16] E. T. Hale, W. Yin, and Y. Zhang. Fixed-point continuation for
l1-minimization: Methodology and convergence. SIAM Journal on
Optimization, 19(3):1107–1130, 2008.
[17] H. W. Haussecker and D. J. Fleet. Computing optical flow withphysical models of brightness variation. IEEE Transactions on
Pattern Analysis and Machine Intelligence (PAMI), 23(6):661–673,
2001.
[18] B. K. P. Horn and B. G. Schunck. Determining optical flow. Artif.
Intell., 17(1-3):185–203, 1981.
[19] C. Lei and Y.-H. Yang. Optical flow estimation on coarse-to-fineregion-trees using discrete optimization. In ICCV, 2009.
[20] V. Lempitsky, S. Roth, and C. Rother. Fusionflow: Discretecontinuous optimization for optical flow estimation. In CVPR, 2008.
[21] C. Liu, J. Yuen, A. Torralba, J. Sivic, and W. T. Freeman. Sift flow:
Dense correspondence across different scenes. In ECCV, volume 3,
pages 28–42, 2008.
[22] D. G. Lowe. Distinctive image features from scale-invariant keypoints. International Journal on Computer Vision (IJCV), 60(2):91–
110, 2004.
[23] B. D. Lucas and T. Kanade. An iterative image registration
technique with an application to stereo vision. In International Joint
Conference on Artificial Intelligence (IJCAI), pages 674–679, 1981.
[24] É. Mémin and P. Pérez. Hierarchical estimation and segmentationof dense motion fields. International Journal on Computer Vision(IJCV), 46(2):129–155, 2002.
[25] S. M.Seitz and S. Baker. Filter flow. In ICCV, 2009.
[26] S. Roth and M. J. Black. On the spatial statistics of optical flow. In
ICCV, pages 42–49, 2005.
[27] C. Rother, V. Kolmogorov, V. S. Lempitsky, and M. Szummer.
Optimizing binary mrfs via extended roof duality. In CVPR, 2007.
[28] P. Sand and S. J. Teller. Particle video: Long-range motion estimation
using point trajectories. In CVPR, volume 2, pages 2195–2202, 2006.
[29] D. Scharstein and C. Pal. Learning conditional random fields forstereo. In CVPR, 2007.
[30] L. Sigal, A. O. Balan, and M. J. Black. Humaneva: Synchronizedvideo and motion capture dataset and baseline algorithm for evaluation of articulated human motion. International Journal on Computer
Vision (IJCV), 87(1-2):4–27, 2010.
[31] M. Sizintsev and R. P. Wildes. Efficient stereo with accurate 3-dboundaries. In BMVC, pages 237–246, 1996.
[32] F. Steinbrücker and T. Pock. Large displacement optical flow
computation without warping. In ICCV, 2009.
[33] D. Sun, S. Roth, and M. J. Black. Secrets of optical flow estimationand their principles. In CVPR, 2010.
[34] D. Sun, S. Roth, J. P. Lewis, and M. J. Black. Learning optical flow.
In ECCV, volume 3, pages 83–97, 2008.
[35] Y. Wang, J. Yang, W. Yin, and Y. Zhang. A new alternating
minimization algorithm for total variation image reconstruction.
SIAM Journal on Imaging Sciences, 1(3):248–272, 2008.
[36] A. Wedel, D. Cremers, T. Pock, and H. Bischof. Structure- andmotion-adaptive regularization for high accuracy optic flow. In
ICCV, 2009.
[37] A. Wedel, T. Pock, C. Zach, H. Bischof, and D. Cremers. An
improved algorithm for tv-l1 optical flow computation. In Dagstuhl
Visual Motion Analysis Workshop, 2008.
[38] M. Werlberger, T. Pock, and H. Bischof. Motion estimation with
non-local total variation regularization. In CVPR, 2010.
[39] M. Werlberger, W. Trobin, T. Pock, A. Wedel, D. Cremers, and
H. Bischof. Anisotropic huber-l1 optical flow. In BMVC, 2009.
[40] J. Xiao, H. Cheng, H. S. Sawhney, C. Rao, and M. A. Isnardi. Bilateral filtering-based optical flow estimation with occlusion detection.
In ECCV, volume 1, pages 211–224, 2006.
[41] L. Xu, J. Chen, and J. Jia. A segmentation based variational model
for accurate optical flow estimation. In ECCV, volume 1, pages
671–684, 2008.
[42] L. Xu, J. Jia, and Y. Matsushita. Motion Detail Preserving Optical
Flow Estimation. In CVPR, 2010.

14[43] C. Zach, T. Pock, and H. Bischof. A duality based approach forrealtime tv-l1 optical flow. Pattern Recognition (Proc. DAGM), pages
214–223, 2007.
[44] H. Zimmer, A. Bruhn, and J. Weickert. Optic flow in harmony.
International Journal on Computer Vision (IJCV), 93(3):368–388,
2011.
[45] H. Zimmer, A. Bruhn, J. Weickert, B. R. Levi Valgaerts and,
Agustı́n Salgado, and H.-P. Seidel. Complementary optic flow.
In International Conference on Energy Minimization Methods in
Computer Vision and Pattern Recognition (EMMCVPR), 2009.

Li Xu received the BS and MS degrees
in computer science and engineering from
Shanghai JiaoTong University in 2004 and2007 respectively, and the PhD degree in2010 in computer science and engineeringfrom the Chinese University of Hong Kong,
where he is currently a postdoctoral fellow.
He received the Microsoft Research Asia
Fellowship Award in 2008 and served as
reviewers for several major computer visionand graphics conferences and journals. His
research interests include motion estimation, motion deblurring, image/video analysis and enhancement. He is a member of the IEEE.

Jiaya Jia received the PhD degree in Computer Science from Hong Kong University of
Science and Technology in 2004 and is currently an associate professor in Departmentof Computer Science and Engineering at the
Chinese University of Hong Kong (CUHK).
He was a visiting scholar at Microsoft Research Asia from March 2004 to August
2005 and conducted collaborative research
at Adobe Systems in 2007. He leads theresearch group in CUHK, focusing specifically on computational photography, 3D reconstruction, practicaloptimization, and motion estimation. He serves as an associate
editor for TPAMI and as an area chair for ICCV 2011. He was on the
program committees of several major conferences, including ICCV,ECCV, and CVPR, and co-chaired the Workshop on Interactive
Computer Vision, in conjunction with ICCV 2007. He received theYoung Researcher Award 2008 and Research Excellence Award
2009 from CUHK. He is a senior member of the IEEE.

Yasuyuki Matsushita received his B.S.,
M.S. and Ph.D. degrees in EECS from the
University of Tokyo in 1998, 2000, and 2003,
respectively. He joined Microsoft Research
Asia in April 2003. He is a Lead Researcher
in Visual Computing Group of MSRA. His
major areas of research are computer vision(photometric techniques, such as radiometriccalibration, photometric stereo, shape-fromshading), computer graphics (image relighting, video analysis and synthesis). Dr. Matsushita served as an Area Chair for IEEE Computer Vision and
Pattern Recognition (CVPR) 2009 and International Conference on
Computer Vision (ICCV) 2009, and he is on the editorial board
member of International Journal of Computer Vision (IJCV), IPSJ
Journal of Computer Vision and Applications (CVA), The VisualComputer Journal, and Encyclopedia of Computer Vision. He serves
as a Program Co-Chair of PSIVT 2010, Demo Co-Chair for ICCV2011, and Program Co-Chair for 3DIMPVT 2011. He is appointed asa Guest Associate Professor at Osaka University (April 2010-) andNational Institute of Informatics, Japan (April 2011-). He is a senior
member of IEEE.