3970 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 29, 2020
Abstract—"is paper proposes graph Laplacian regularization
for robust estimation of optical flow. First, we analyze the spectral
properties of dense graph Laplacians and show that dense graphs
achieve a better trade-off between preserving flow discontinuities
and filtering noise, compared with the usual Laplacian. Using this
analysis, we then propose a robust optical flow estimation method
based on Gaussian graph Laplacians. We revisit the framework of
iteratively reweighted least-squares from the perspective of graph
edge reweighting, and employ the Welsch loss function to preserve
flow discontinuities and handle occlusions. Our experiments using
the Middlebury and MPI-Sintel optical flow datasets demonstrate
the robustness and the efficiency of our proposed approach.
Index Terms—Optical flow, regularization, optimization, robust
estimation.
I. I!"#$%&'"($!
N IMAGE PROCESSING and machine vision, optical flow
refers to a pattern of visible motion of objects, surfaces, and
edges in a visual scene, created by the relative motion between
the observer and the scene. Optical flow is often estimated as a
first step in a number of imaging- and vision-related tasks such
as image registration, object detection and tracking, navigation
and visual odometry for robots.
Posed as a mathematical optimization problem, optical flow
estimation exhibits certain similarities to the problem of image
denoising, deblurring and inpainting, which are all examples of
the discrete inverse problems typically encountered in the field
of image processing. At the same time, however, the estimation
of optical flow differs from the two related inverse problems in
two major aspects. First, optical flow estimation is non-convex
due to the nonlinearity in the image itself, and this differs from
the image deblurring or the inpainting problems, both of which
can be formulated as convex problems. Second, the estimation
involves using spatiotemporal derivative images which contain
large amounts of noise.
Besides the variational approaches [1], [2], feature-matching
based approaches [3], [4] have been proposed more recently to
address complications due to image non-linearity. Probabilistic
frameworks based on message passing [5], [6] aim similarly to
resolve this non-convexity, but by discretization of the original
continuous optimization problem, and optimizing directly over
the space of discrete flow vector labels. In this paper, we focus
our discussion on variational approaches, and, in particular, on
those resembling the original one due to Horn and Schunck [7]
since our method can be seen as an extension of [7] to a graph-
like space. Fundamentally, all variational approaches minimize
an energy function consisting of data fidelity and regularization
terms, using the machinery of numerical methods to efficiently
find the solution. Methods for estimating larger motion rely on
feature matching as a first step [3], [6], [8] and perform optical
flow estimation to further refine the large motion. Zis suggests
that optical flow estimation is indispensable to estimating large
motion as well. Recent self-supervised learning methods based
on deep convolutional neural networks [9], [10] can learn more
complex variational models directly from image sequences.
Some time ago, Gilboa and Osher [11] examined Tikhonov-
regularization based on nonlocal operators with applications to
inverse problems in vision. Zey define the so-called non-local
gradient, divergence and Laplacian operators in the continuous
domain, analogously with their graph-based counterparts. Zeir
original justification for the nonlocal operators in [11] involves
the problems of image texture denoising and inpainting, where
the imposition of regularity on texture (an inherently non-local
property) makes sense.
Interestingly, non-local operators have also been found to be
useful for the reconstruction of piecewise smooth signals, such
as optical flow fields, with the non-local operators appearing in
the form of weighted median filtering [1], a bilateral filter term
[12] and nonlocal total variation regularization [2], [13]. All of
them derive soft segmentation cues from the color pixel values
to determine the support of the nonlocal kernels. Ze empirical
results from [1], [2], [13] show the advantages of the non-local
gradient, divergence and Laplacian operators over the classical
discretized counterparts. However, texture is not characteristic
of piece-wise smooth signals, and a different justification needs
to be given for the advantages of nonlocal regularization in the
reconstruction of piece-wise signals.
In this work, we show that the benefit of non-local operators
in recovering piecewise smooth signals comes mostly from the
spectra of their corresponding dense graph Laplacians. We first
examine the properties of the Laplacian matrices of such dense
graphs to reveal a number of important results concerning their
graph spectra and Green’s functions. Based on the analysis, we
then propose our robust optical flow estimation method using a
Gaussian (bilateral) graph Laplacian regularizer, also revisiting
Graph Laplacian Regularization for
Robust Optical Flow Estimation
Sean I. Young , Member, IEEE, Aous T. Naman , Member, IEEE, and David Taubman , Fellow, IEEE
I
)*+,-./0123 /4.405463 %4.4784/3 9:3 ;<=>?3 /450-463 @*+,*/A3 B:3 ;<=B3 *+633
)*A3;B:3;<=B?3*..412463C4124784/3;D:3;<=BE3F43*--G.0*2434602G/3.GG/60+*20+H3
2I43 /4504J3 GK3 2I0-3 7*+,-./0123 *+63 *11/G50+H3 023 KG/3 1,8L0.*20G+3 J*-3 M/GKE33
'I,+70+H3N0E3!"#$$%&'#()*(+,-./0#$1,2%-(,34,5#.(+46,
CE3(E3OG,+H3J*-3J02I32I43C.IGGL3GK3PL4.2/0.*L3P+H0+44/0+H3*+63"4L4.G77,Q
+0.*20G+-:3F43&+054/-02A3GK3!4J3CG,2I3R*L4-:3CA6+4A:3!CR3;<S;3T,-2/*L0*E3
U430-3+GJ3J02I32I43%41*/274+23GK3PL4.2/0.*L3P+H0+44/0+H:3C2*+KG/63&+054/-02A:3
C2*+KG/6:3'T3BV9<S3&CT3W4Q7*0LX3-4*+<Y-2*+KG/6E46,ZE3
TE3"E3!*7*+3*+63%E3"*,87*+3*/43J02I32I43C.IGGL3GK3PL4.2/0.*L3P+H0+44/0+H3
*+63"4L4.G77,+0.*20G+-:3F43&+054/-02A3GK3!4J3CG,2I3R*L4-:3CA6+4A:3!CR3
;<S;:3T,-2/*L0*3W47*0LX3*G,-Y,+-JE46,E*,?36E2*,87*+Y,+-JE46,E*,ZE3
%0H02*L3$8[4.23(64+20\4/3=<E==<B]"(ME;<=BE;BVSDS93
=<S>Q>=VB3^3;<=B3(PPPE3M4/-G+*L3,-430-314/702246:38,23/41,8L0.*20G+]/460-2/08,20G+3/4_,0/4-3(PPP314/70--0G+E3
C443I221X]]JJJE0444EG/H]1,8L0.*20G+-`-2*+6*/6-]1,8L0.*20G+-]/0HI2-]0+64aEI27L3KG/37G/430+KG/7*20G+E3
YOUNG et al.: GRAPH LAPLACIAN REGULARIZATION FOR ROBUST OPTICAL FLOW ESTIMATION 3971
the framework of iteratively reweighted least-squares from the
perspective of graph edge reweighting. Ze robust, Welsch loss
function is incorporated to preserve flow discontinuities, and to
handle occlusions. Since the Welsch loss is nonconvex, we also
propose a graduated nonconvexity scheme for our Welsch loss-
based optimization objective.
Ze reader may wish to refer to [14], [15] for an overview of
graph spectral filtering and its relationship to regularization.
II. OM"('TN3FN$R3Eb&T"($!C
Var i at i ona l o pt ic a l fl ow est ima tio n me th ods s eek to c om pu te
the per-pixel motion across two picture frames captured at two
different time instances. Let us denote by
0
and
1
∈
𝑁
, the
raster ordering of the two discrete picture frames. Ze so-called
brightness consistency equations, written in matrix form are
𝑥
𝑥
+
𝑦
𝑦
+
𝑡
= ,
(1)
in which
∗
= (
∗
), and
𝑥,𝑦,𝑡
denote the spatiotemporal
derivatives of picture
0
in the -, - and -directions. Ze two
flow vector components 
𝑥
,
𝑦
cannot be computed without
introducing additional assumptions since (1) is a rank-deficient
problem as is seen from the fact that (1) has double the number
of unknowns, the elements of 
𝑥
,
𝑦
, as there are equations.
Ze rank-deficiency of problem (1), however, is not the only
source of “ill-conditioning”. Suppose the flow was constrained
to lie along the -direction, as is usually the case with disparity
estimation. In this case, we have the system of equations
𝑥
𝑥
+
𝑡
= ,
𝑦
= ,
(2)
one that is equally ill-posed, even with the “correct” number of
equations in the unknowns—an arbitrary row of the coefficient
matrix
𝑥
can happen to be zero to render problem (2) devoid
of a unique solution. Furthermore, the error introduced into the
coefficients
𝑥
from numerical differentiation can amplify the
high-frequency noise in the naïve solution of (2), often referred
to as the inversion noise in discrete inverse problems [16].
In either case, the ill-conditioning requires us to reformulate
or regularize the problem for numerical treatment. To facilitate
the presentation, we use problem (2) as the running example in
our analysis. Ze graph-based regularization which we propose
later is rotationally invariant and separable across the two flow
components, and this allows us to limit the analysis to only one
component. We consider the solution to problem (1) during the
actual numerical optimization (Section X).
III. OM"('TN3FN$R3EC"()T"($!
Consider the non-parametric formulation of the optical flow
estimation problem, in which we estimate a flow vector at each
pixel location subject to a quadratic regularity constraint:
(3)
in which is a Laplacian matrix (again, we consider the stereo
case first for analysis). Ze constrained problem (3) is typically
relaxed using the method of Lagrange multipliers. Ze Karush-
Kuhn-Tucker conditions guarantee that the solution of problem
(3) is obtained alternatively by solving the relaxation:
minimize () =
‖
𝑥
+
𝑡
‖
2
2
+
2
∗
,
(4)
in which
2
is our dual variable for the inequality constraint in
problem (3). Ze relaxation (4) is essentially the formulation of
Horn and Schunck [7]. It is known that (4) tends to smooth the
flow discontinuities excessively due to the isotropy of .
To better investigate the role of , it is instructive to express
the data term of (4) in the normalized form [17] to yield
minimize () =
‖
+
𝑥
†
𝑡
‖
2
2
+
2
∗
,
(5)
to write problem (4) as one of denoising the transformed signal
= −
𝑥
†
𝑡
, where
†
is the pseudo-inverse of . Since is
the quotient of two noisy high-pass signals
𝑡
and
𝑥
, one may
expect to be very noisy in general, with a low signal-to-noise
ratio. Ze quality of the solution of problem (5) therefore relies
upon the spectral properties of the Laplacian . If we appeal to
the eigen-decomposition = 
∗
of the Laplacian, we can
express the solution of (5) as
opt
= (+
2
)
−1
∗
,
(6)
that is, spectral filtering of using the spectral filter factors =
(
1
,
2
, . . . ,
𝑁
), where
𝑛
=
1
1 +
2
𝑛
, for 1 ≤≤,
(7)
so the filter factor for the th graph frequency depends equally
on
𝑛
and the global smoothness parameter
2
.
Zerefore, our goal is now to construct a Laplacian operator
(a) (b) (c) (d)
c0HE3=E3d/*1I-3GK360e4/4+2364+-02A:30+6,.4638A360e4/4+232A14-3GK3N*1L*.0*+3f4/+4L-3-IGJ+30+38GL6E3W*Z3-4.G+6Q60e4/4+.43f4/+4L:3W8Z340HI2Q1G0+23N*1L*.0*+3f4/+4L:3
W.Z3+G+LG.*L3;VQ1G0+23N*1L*.0*+3f4/+4L:3*+63W6Z3+G+LG.*L3gGJQH,06463G/307*H4QH,06463;VQ1G0+23N*1L*.0*+3f4/+4LE3
3972 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 29, 2020
that provides an improved trade-off between filtering noise and
preserving flow discontinuities. We see later that the Laplacian
of a dense graph provides the desired trade-off behavior. In our
paper, a graph refers to a weighted graph, and is denoted using
= (, , ) with vertices , edges and a weight function
∶→
+
. Given some graph , we refer to its vertices and
edges using () and () respectively.
IV. E(dP!hTN&PC3$c3DP!CP3LTMNT'(T!C
A. Laplacians vs Graph Laplacians
Ze discrete Laplacian operator , typically used in classical
Tikhonov-regularized inverse problems can be expressed using
a filter kernel. Zis kernel is
(8)
that is, a linear operator which takes the difference between its
input and a low-pass filtered version of it. Ze two-dimensional
analog of (8) is considered in [7]. Since the space
𝑁
⊂
+
𝑁
of
Laplacian matrices has a bijection with the space
𝑁
of graphs
on vertices, we can associate (8) with the grid graph
𝑟=1
shown in Fig. 1 (a). Ze two-dimensional version
of this graph
is also depicted in Fig. 1 (b). From a graph perspective, we can
understand the robust optical flow estimation method of [18] as
locally modifying the isotropic kernel (8) to
2
R
(
)
=
[
0
2
0
]
−
2
2 + 2
[
1 1 +
]
(9
)
(where 0 ≤≤1) or
2
L
() =
[
0
2
0
]
−
2
2 + 2
[
+ 1 1
]
(10)
depending on whether a flow discontinuity was detected on the
right- or the left-hand side at the given location in the previous
flow estimate.
Later, we see that if contains both flow discontinuities and
noise, localized Laplacian kernels of the form (8)–(10) provide
poor trade-offs between preserving the flow discontinuities and
filtering noise in signal —using a smaller value of provides
inadequate filtering of the noise while a larger value of blurs
the flow discontinuities. Since the optical flow we can estimate
is largely translational (i.e. constant) at least locally, it is better
to filter out most frequencies except around DC, where smooth
components of the optical flow field belong. However,
𝑟=1
is
poor at filtering lower-frequencies (Fig. 2, left)—one could use
a large global value of (7) to boost the
𝑛
, but this blurs flow
discontinuities unless = 0 in all instances of (9) and (10). On
the other hand, denser graphs
𝑟=4
and
𝑟=16
produce a more
uniform set of eigenvalues, so one can use a small global value
of for the same amount of denoising, while better preserving
the flow discontinuities.
Let us define
𝑟
for = 1, 2, . . ., as weighted graphs whose
graph vertices correspond to the locations
𝑛
=
𝑛
𝑛
, 1 ≤≤,
(11)
and the set of edges given by
(
𝑟
) = {(
𝑛
,
𝑚
) ∈
2
∶
|
𝑛
−
𝑚
|
∞
≤1},
(12)
in which | · | denotes a norm on
2
. A graph
𝑟
with a larger
contains more edges for the same number of vertices. We refer
to such a graph as being dense, and illustrate it in Fig. 1 (c) and
(d). Ze Laplacian of
𝑟
can be defined as = −, where
and are respectively the degree and adjacency matrices of
the graph
𝑟
.
B. Box-weighted Graph Laplacians
We also need to define a weight function for the edges of
𝑟
in terms of distance between adjacent vertices. In our work, we
consider the box and the Gaussian weights, both of which have
theoretical and practical significance. Ze two weight functions
perform similarly in practice, but the Gaussian weight function
= 1
= 2
= 4
= 16
= 1
= 2
= 8
= 0.5
c0HE3;E3P0H4+5*L,4-3𝜆
𝑚
:3𝑚 = 0,1, . . . ,400:3GK32I43N*1L*.0*+3𝐋 ∈ 𝕊
+
800
31LG22463KG/32I438GaQJ40HI2463WL4K2Z3*+63d*,--0*+QJ40HI2463W/0HI2Z3.*-4-E3F438GaQJ40HI2463
H/*1I3N*1L*.0*+-31/G6,.43G-.0LL*20+H340H4+5*L,4-3JI0L432I4340H4+5*L,4-3𝜆
𝑚
3GK32I43d*,--0*+QJ40HI2463H/*1I3N*1L*.0*+-364./4*-437G+G2G+0.*LLA3*-3
𝑚 → 0E3"JGQ
6074+-0G+*L3-14.2/*3-IGJ+30+32I430+-42-3.G//4-1G+632G3𝑟 = 2
3W8GaZ3*+63
𝜎 = 1
3Wd*,--0*+ZE3F438GaQJ40HI246340H4+5*L,4-3*/43+G23.0/.,L*/LA3-A7742/0.E
YOUNG et al.: GRAPH LAPLACIAN REGULARIZATION FOR ROBUST OPTICAL FLOW ESTIMATION 3973
renders a Gaussian filter, a property that can be exploited for
efficient evaluation of the mapping . To show the effect
of graph density better, all our graphs are degree-normalized to
unity. First, we analyze the Laplacians of graphs weighted with
the box function
𝑛𝑚
=
1 if
𝑛
∼
𝑚
,
0 otherwise,
(13)
corresponding also to the case of an un-weighted graph. Dense
graphs produce more uniform eigenvalues as can be seen from
the following proposition.
Proposition 1. Ze eigenvalue
𝑚
for the th Fourier mode or
the eigenvector of the graph
𝑟
on vertices is given by
𝑚
=
−1
1 −
sinc(
⁄
)
sinc
(
⁄
)
,
(14)
for all = 1, 2, . . . , , where = 2+ 1 is the degree of the
graph
𝑟
.
Proof. Ze standard results from spectral graph theory provide
the bounds on
𝑚
[19]. A general expression for the spectrum
of circulant graphs is derived in [20]. One can prove (14) using
the fact that the eigenvectors of a circulant
3
∈
𝑁×𝑁
are the
discrete Fourier vectors. Zerefore, the eigenvalues of the graph
𝑟
are given by the -point DFT of the sequence
𝑟
[
]
= (1 + 1 2
⁄
)
[
]
−(1 2
⁄
)Π
𝑁,2𝑟+1
[
]
,
(15)
in which Π
𝑁,𝑙
denotes the -periodic box sequence of length
symmetric about = 0. Ze magnitude responses of Π
𝑁,𝑑
and
in (15) can be derived using the standard results from digital
signal processing [21]. Linearly combining the two magnitude
responses as per (15) produces (14).
In particular, one can reproduce the eigenvalues of the usual
second difference operator up to scale, by substituting = 3 in
(15). As = 2+ 1 →, the eigenvalues
𝑚
of
𝑟
approach
1
= 0,
𝑚
=
−1
for all = 2, 3, . . . ,
(16)
known as the spectrum of the “complete graph”
𝑁
in spectral
graph theory [19, Sec. 1.2].
Ze spectra of box-weighted graphs are illustrated in the left-
hand side of Fig. 2 for varying graph densities. We see that the
spectra of graphs of increasing density approach the brick-wall
spectrum of the complete graph. Zis suggests that too dense a
graph can unduly penalize non-constant (but smooth) trends in
the flow.
C. Gaussian-weighted Graph Laplacians
Our treatment of weighted graph Laplacians can be extended
to the Gaussian weighted case. Since it is more natural to use a
scale parameter to parameterize a Gaussian, we now denote by
𝜎
, a weighted graph whose vertices
𝑛
∈(
𝜎
) correspond
to the locations
𝑛
=
𝑛
𝑛
, 1 ≤≤,
(17)
the set of edges given by
(
𝜎
)
= {
(
𝑛
,
𝑚
)
∈
2
∶
|
𝑛
−
𝑚
|
2
≤2},
(18)
and the weighting function takes the form
(
𝑛
,
𝑚
) = e
−|𝐩
𝑛
−𝐩
𝑚
|
2
/2
.
(19)
Ze Gaussian weighting function (19) is seen in applications
such as bilateral filtering [22], non-local means denoising [23]
and kernel density estimation. Ze fast Gauss transform used in
these applications can similarly be exploited in our problem to
estimate the optical flow efficiently. We later appeal also to the
relationship between the Welsch penalty [24] and the Gaussian
function for a robust formulation of the optical flow estimation
problem.
In this analysis, we assume the Gaussian weights are applied
to every edge of a complete graph. In practice, edges that have
insignificant weights can be pruned as in (18). Ze eigenvalues
of the Gaussian-weighted graph Laplacian can be given via the
following proposition.
Proposition 2. Ze eigenvalue
𝑚
for the th Fourier mode or
the eigenvector of the graph
𝜎
on vertices is given by
𝑚
=
−1
−
3
2
⁄
, e
−1 2𝜎
2
⁄
,
(20)
for all = 1, 2, . . . , , where
3
is the Jacobi theta function
defined as
3
(, ) = ∑
𝑡=1
∞
2
𝑡
2
cos() + 1
(21)
[25, Sec. 16.27], and is the degree of
𝜎
given by
=
3
0, e
−1 2𝜎
2
⁄
.
(22)
Proof. Ze rows in the Gaussian weighted Laplacian matrix are
produced by the -periodic sequence
𝜎
[
]
= (
[
]
−exp(−
2
2
2
⁄
)) (−1)
⁄
(23)
in which
= ∑
𝑡=−∞
∞
exp
(
−
2
2
2
⁄
)
(24)
is the degree of
𝜎
. Observe that lim
𝜎→∞
/=
√
2 and the
frequency response of
𝜎
is the sum of those of the Dirac delta
and the sampled Gaussian sequence.
Ze standard results from Fourier transforms [26] tell us that
the frequency response of a sampled signal is given by periodic
summation. Ze response of a sampled Gaussian is therefore
[
]
= ∑
𝑡=−∞
∞
exp(−2
2
2
(
⁄
−)
2
),
(25)
and using the definition of the Jacobi theta function
3
to write
both (24) and (25) produces the required result.
Ze eigenvalues of
𝜎
are graphed on the right-hand side of
Fig. 2 for different values of . In contrast with the eigenvalues
of dense box-weighted graph Laplacians where Gibbs’ effect is
apparent, the eigenvalues of the Gaussian-weighted Laplacians
decrease monotonically as →0, which is expected from the
Gaussian function itself. Zis property is optimal since smooth
modes are never regularized more strongly than the oscillatory
ones. Ze circular symmetry and separability of Gaussians also
3974 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 29, 2020
lead to efficient evaluation of the mapping → with a fast
Gaussian filter implementation.
V. TUP3G#PP!3c&!'"($!C
To reduce the amount of blur across flow discontinuities, we
additionally require a mechanism to de-weight the graph edges
that straddle potential flow interfaces. In flow estimation, such
de-weighting can be achieved using image-guided [27], [28] or
flow-guided [29]–[31] approaches, which change edge weights
based on discontinuities in the image texture, or the previously
estimated flow field, respectively. We show one example of the
resulting graph in Fig. 1 (d). In the remainder of this paper, we
focus on Gaussian-weighted Laplacians parameterized by .
To observe the behavior of denser graph Laplacians across a
flow discontinuity, consider problem (5) again, noting this time
its solution
opt
is linear in = −
𝑥
†
𝑡
. Said differently,
opt
= , = (+
2
)
−1
,
(26)
and such a matrix is referred to as the Green’s function [32]
of problem (5). In spectral graph theory, is referred to as the
doubly-stochastic graph matrix [33] of .
In particular, one can regard the th row of as prescribing
a particular convex combination of the elements of the signal
to produce the th element of
opt
. We also see from (6) that a
circulant Laplacian produces a circulant , so we can obtain
the solution
opt
by convolving with a column of . We plot
in Fig. 3, the convolution kernels of graph matrices of different
density.
If certain graph edges are de-weighted, the resulting Green’s
function is no longer circulant nor diagonalized by the discrete
Fourier vectors. We derive a general expression for the Green’s
function of (5) in terms of a circulant matrix and successive
rank-1 updates to it.
Proposition 3. Given a circulant graph on vertices which
has the weight function , let us define a new graph
, having
the weight function
(
𝑛
,
𝑚
) =
1(
𝑛
,
𝑚
) if
𝑛
∼
𝑈
𝑚
,
(
𝑛
,
𝑚
) otherwise,
(27)
in which
{
,
𝑐
}
is some non-trivial partition of (). Let us
also denote by ∈
𝑁×𝑀
, the (weighted) incidence matrix of
the subgraph of
whose edges have been de-weighted. We
can express the Green’s function
𝑀
of
recursively as
𝑚
=
𝑚−1
+ ∆
𝑚−1
, 0 < ≤,
(28)
and
0
is the Green’s function of , while
∆
𝑚−1
=
(
𝑚−1
𝑚
)(
𝑚−1
𝑚
)
∗
(
1 −
)
−1
−2
−
𝑚
∗
𝑚−1
𝑚
,
(29)
and
𝑚
is the th column of .
Proof. Observe that we can write

∗
=
1
1
∗
+
2
2
∗
++
𝑀
𝑀
∗
,
(30
)
from which, a successive application of the Sherman-Morrison
formula [34] yields a straight-forward proof by induction.
In particular, one can interpret the th row of 
𝑚−1
as an
adjustment of the coefficients producing the th element of the
solution
opt
of problem (5). Zis adjustment is proportional to
(
𝑚−1
𝑚
)
∗
, the slope
1
of the rows of the unadjusted Green’s
function
𝑚−1
. Since each
𝑚−1
is doubly-stochastic, we see
that 
𝑚−1
has a zero row-sum. Zerefore, 
𝑚−1
transfers
part of the weight mass of
𝑚−1
from the incorrect side of the
partition to the correct side. Such an adjustment increases with
the number of edges that cross the discontinuity, so a dense
graph effects a larger and hence better adjustment than a sparse
one. Fig. 4 shows that for a similar amount of denoising, dense
graphs smooth less across soft discontinuities, which preserves
flow discontinuities better.
1
C0+.43𝐡
𝑚
0-3*3-1*/-4354.2G/3.G+2*0+0+H3−
√
𝑤3*+63
√
𝑤3J02I3-G743ℎ3i4/G-30+Q
842J44+:3023LG.*LLA360e4/4+20*24-32I43/GJ-3GK3𝐆
𝑚−1
3*23-.*L43ℎE3
=
35
=
60
=
120
=
4
=
8
=
16
c0HE39E3F434002I3/GJ3𝐠3GK32I4360-./4243d/44+3K,+.20G+3𝐆3GK31/G8L473W;DZ30+32I43G+4Q6074+-0G+*L3.*-43JI4+3𝐺
𝜎
30-3/4H,L*/E3MLG2-3H054+3.G//4-1G+632G32I43.*-4-3GK3
-1*/-43𝐺
𝜎=0.5
3WL4K2Z3*+63
𝐺
𝜎=8
3W/0HI2Z3H/*1I-E3F43/GJ-3GK3
𝐆
3*/43NG/4+2i0*+QL0f43K,+.20G+-E
YOUNG et al.: GRAPH LAPLACIAN REGULARIZATION FOR ROBUST OPTICAL FLOW ESTIMATION 3975
VI. RPRP(dU"P%3LTMNT'(T!C
A majority of flow discontinuities are discovered only as the
estimation progresses. Zis is the case e.g. at the boundaries of
two similarly colored objects under different motion. Since the
graph Laplacian applies strong quadratic regularization over
all graph edges, we require a procedure which can reweight the
graph edges iteratively, based on the previous flow estimate. In
the PDE-based imaging formalism, this procedure is known as
the lagged diffusivity method [35]. In robust estimation, this is
the method of iteratively reweighted least-squares [36].
A. Reweighted Laplacians for Robustness
To relate edge reweighting to robust flow estimation, we can
generalize the quadratic formulation (5) to
minimize () =
‖
+
𝑥
†
𝑡
‖
2
2
+
2
‖
(
∗
)
‖
1
1
,
(31)
denoting by , the incidence matrix of . With () =
2
, we
recover our original problem (5). Problem (31) becomes robust
if the function ∶→
+
increases sub-quadratically, i.e., the
derivative ′ is bounded. Ze loss functions we consider can be
classified as either Huber-type (Huber, and truncated-quadratic
losses) or Charbonnier-type (Charbonnier, Lorentzian, Geman-
McClure and Welsch losses). Zeir exact definitions [37, p. 24]
are unimportant for our discussions.
Associated with any loss function is its weighting function
() = ′()/. We can solve the non-quadratic problem (31)
using the fixed-point iteration
𝑗+1
= −+
2

𝑗
∗
−1
𝑥
†
𝑡
(32)
for = 0,1, . . . in which
𝑗
= ((
∗
𝑗
)) are the lagged
diffusivity weights obtained using the th flow estimate
𝑗
. We
can equivalently view
𝑗+1
as the minimizer of the graph edge
reweighted least-squares problem
minimize
𝑗+1
() =
‖
+
𝑥
†
𝑡
‖
2
2
+
2
∗
𝑗
,
(33)
in which we denote
𝑗
= 
𝑗
∗
for short.
Ze equivalence between (32) and solving (33) may be used
to combine the reweights arising from the iterative process and
the static graph edge weights. Suppose that the weight function
() = exp(−
2
2
2
⁄
), and our Laplacian is also weighted in
the spatial domain with a Gaussian (19). Zen,
𝑗
corresponds
to the Laplacian of a graph whose vertices are
𝑛
=
𝑛
𝑋
𝑛
𝑌
𝑛
𝑈
, 1 ≤≤,
(34)
in which
𝑛
is the th element of
𝑗
, and the edge weighting is
given by (
𝑛
,
𝑚
) = e
−|𝐩
𝑛
−𝐩
𝑚
|
2
/2
as before.
B. Gauss Meets Welsch
Ze Gaussian is in fact the weighting function of the Welsch
loss. It is also the limit of the weight functions of Charbonnier-
type losses as the degree of non-convexity of the loss functions
increases. To see this, we generalize the weighting functions of
Charbonnier-type losses [37, p. 24] as
(𝜂)
() = (1 +
−1
2
2
2
⁄
)
−𝜂
,
(35)
the choices = 0.5, 1 and 2 corresponding, respectively, to the
Charbonnier, Lorentzian and Geman-McClure losses. We have
the weighting function of the Welsch loss (the Gaussian) in the
limit of
(𝜂)
as →∞ (this can be proved using the definition
of Euler’s number, e). By contrast, the weighting functions for
the Huber-type losses converge to a box function. Both Huber-
and Charbonnier-type weighting functions are shown in Fig. 5.
In practice, the convex Huber or Charbonnier loss functions
corresponding to the decay rate = 0.5 are inadequate to curb
the influence of gross outliers [36], which, for our optical flow
estimation problem, correspond to flow discontinuities or large
data term errors in the occluded regions, suggesting that a non-
convex loss function be used instead. However, the Huber-type
loss functions may lead to a non-smooth optimization behavior
=
35
=
60
=
120
=
4
=
8
=
16
c0HE3VE3F434002I3/GJ3GK360-./4243d/44+3K,+.20G+3𝐆3GK31/G8L473W;DZ30+32I43G+4Q6074+-0G+*L3.*-43JI4+3𝐺
𝜎
30-30//4H,L*/E3RI4+3J43I*543*3-GK2360-.G+20+,02A3842J44+3
2I43-,8-42-3𝑈 = {1,2, . . . ,400}
3*+63
𝑈
𝑐
3GK354/20.4-3
𝑉 ( 𝐺
𝜎
) = {1,2, . . . ,800}
:32I43-1*/-43H/*1I3
𝐺
𝜎=0.5
3WL4K2Z3KG/7-32I43*54/*H43*./G--32I4354/20.4-30+38G2I3
𝑈
3*+63
𝑈
𝑐
33
*+63I4+.438L,/-32I4360-.G+20+,02AE3F4364+-43H/*1I3𝐺
𝜎=8
3W/0HI2Z3*54/*H4-37G-2LA3*./G--3G+LA32I4354/20.4-3GK3
𝑈
:384224/31/4-4/50+H3gGJ360-.G+20+,0204-E
3976 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 29, 2020
during reweighted least squares since their weighting functions
are not differentiable. By contrast, Charbonnier-type weighting
functions are infinitely differentiable, which allows the solution
to evolve smoothly. A similar observation is reported in [36].
Associated with Gaussian weighting is the Welsch loss
() =
2
(1 −exp(−
2
2
2
⁄
)),
(36)
the Charbonnier-type equivalent of the truncated quadratic loss
[38]. In this paper, we consider the Welsch loss exclusively for
edge reweighting and occlusion handling, noting that Gaussian
edge reweighting is amenable to a fast filtering implementation
using the fast Gauss transform. To control the non-convexity of
the Welsch function, we additionally propose in Section IX our
graduation scheme inspired by [38, Sec. 7].
VII. G#TMU3LTMNT'(T!3C$!C"#&'"($!
To construct a dense anisotropic graph Laplacian, we appeal
to the concept of Gestalt grouping introduced in [39], also used
by several authors [1], [2] for estimating optical flow. Similarly
colored pixels inside a neighborhood most likely correspond to
a part of the same object undergoing similar motion; also, flow
discontinuities can only be observed across intensity (or color)
change, so it is natural to consider edge weights given in terms
of intensity (and color) difference across adjacent pixels within
a certain locality.
Ze vertices of our graph thus become
𝑛
=
𝑛
𝑋
𝑛
𝑌
𝑛
𝑈
𝑛
𝑍
𝑛
𝑍
𝑛
𝑍
(37)
in which
𝑛
,
𝑛
and
𝑛
represent the CIE Lab values of the th
sample. Our edge weight function (
𝑛
,
𝑚
) = e
−|𝐩
𝑛
−𝐩
𝑚
|
2
/2
as before. Whereas
𝑛
could have involved the RGB values of
the samples instead, we empirically find that the CIE Lab color
space is better-suited for graph construction
2
.
Although our graph is not circulant, we can still interpret the
eigenvectors and eigenvalues of as the graph Fourier vectors
and frequency response [14]. Since is a Gaussianly weighted
graph, the adjacency matrix corresponds to a Gaussian filter
in the ----- space. Noting that = −, the mapping
 can be computed efficiently inside conjugate gradients
using one element-wise multiplication and one application of a
separable, six-dimensional Gaussian filter. In our work, we use
the Gaussian D-tree implementation [40] of the filter.
Ze un-normalized graph Laplacian matrix has a tendency
to regularize a flow vector more strongly at a vertex if many of
its adjacent vertices have similar color or intensity. At the same
time, our would penalize the affine trends in the flow such as
those due to rotation and zoom of the scene. To overcome both
shortcomings, we use the second-order (normalized) Laplacian
operator 
−2
in place of . If is regular, the eigenvalues
of the 
−2
are simply the squares of those of , amplifying
(resp. reducing) the eigenvalues that are greater than (resp. less
than) 1. More importantly, 
−2
is a fourth-order differential
operator so that affine trends in the flow are not penalized. Zis
is equivalently stated as
‖
†

‖
2
2
= 0 for affinely varying .
VIII. HT!%N(!d3O''N&C($!C
In the optical flow problem, the terms “motion” and “optical
flow” are sometimes used interchangeably. Technically, optical
flow refers to motion between a pair of images only in the limit
as the images are sampled infinitesimally close in time with no
occlusions
3
. In practical estimation, however, occluded regions
are likely to be present in the images, due to the low frame rate
of the captured image sequence. Ze need to address occlusion
partly differentiates the optical flow problem from other related
problems such as image deblurring and denoising, in which the
2
(7*H4-3*/432A10.*LLA3-2G/4630+32I43H*77*Q.G//4.2463𝑅′𝐺′𝐵′3KG/7*2:3-G3G+43
7,-23,-432I430+54/-4Q.G//4.2463𝑅𝐺𝐵35*L,4-32G3.G71,2432I43'(P3N*835*L,4-E3
=
0
.
5
=
1
=
2
→∞
=
0
.
5
=
1
=
2
→
∞
c0HE3SE3U,84/Q2A143J40HI23K,+.20G+-3WL4K2Z3.G+54/H432G32I438Ga3K,+.20G+3*-32I4364.*A3/*243GK3J40HI2-3𝜂 → ∞E3F43/*24-3𝜂 = 0.5,1, ∞3.G//4-1G+6:3/4-14.2054LA:32G3
2I43U,84/:3TA5*.03*+632/,+.*2463_,*6/*20.3J40HI2-E3$+32I43G2I4/3I*+6:3'I*/8G++04/Q2A143J40HI23K,+.20G+-3W/0HI2Z3.G+54/H432G32I43d*,--0*+3*-3𝜂 → ∞E3F43/*24-3
𝜂 = 0.5,1,2, ∞
3.G//4-1G+6:3/4-14.2054LA:32G32I43'I*/8G++04/:3NG/4+2i0*+:3d47*+Q).'L,/43*+63d*,--0*+3J40HI2-E
3
"JG3K,/2I4/3*--,7120G+-3*/43*L-G3+44646:3+*74LA3N*784/20*+3/4g4.20G+3*+63
.G+-2*+230LL,70+*20G+:38,232I4-43,-,*LLA3IGL63,+64/370L63.G+6020G+-E
3
YOUNG et al.: GRAPH LAPLACIAN REGULARIZATION FOR ROBUST OPTICAL FLOW ESTIMATION 3977
occlusion phenomenon is not observed.
Larger data term errors encountered in the occlusion regions
can dominate the objective of (5), and we require a mechanism
to deactivate the distortion term in these regions. We can either
assign smaller weights at pixel locations with larger distortions
by iteratively reweighting the distortion term [41], or by jointly
estimating occlusion variables to subsidize larger distortions or
deficits [42]. Both mechanisms have the effect of extrapolating
the unoccluded flow into occlusion regions, as done in [43].
To give an example of the latter approach, Ayvaci et al. [42]
propose to explicitly estimate the occlusion variable by solving
the non-smooth problem
minimize
(
,
)
=
‖
𝑥
+
𝑡
−
‖
2
2
(38)
+
2
∗
+ 2
‖
log(abs() + 1)
‖
1
1
,
using the iteratively re-weighted
1
method [44]. Ze objective
of (38) is to jointly estimate , which subsidizes large intensity
differences in the first term, while enforcing the sparsity of the
subsidizing variable using the last term. However, one could
minimize the objective of (38) over in advance first to obtain
the equivalent problem
minimize () =
‖
(
𝑥
+
𝑡
)
‖
1
1
+
2
∗
,
(39)
in which (see Appendix A for the derivation)
() =
2
if || < 1,
log(+
|
|
+ 1) −log 2
+
(
1 4⁄
)(
−|| −1
)
2
otherwise,
(40)
and =
2
−2
|
|
−3. Ze weighting function of the loss
(40) is plotted as = 1 (up to a scale factor) in Fig. 5 (left).
Problem (39) can now be solved with iteratively-reweighted
least-squares similarly to the edge-reweighted case (31). Let us
denote by
𝑗+1
, the solution of
minimize
𝑗+1
() =
‖
𝑥,𝑗
+
𝑡,𝑗
‖
2
2
+
2
∗
,
(41)
in which we write
𝑥,𝑗
=
𝑗
1 2
⁄
𝑥
and
𝑡,𝑗
=
𝑗
1 2
⁄
𝑡
in terms
of the reweighting matrix
𝑗
= ((
𝑥
𝑗
+
𝑡
)).
As an example of the former approach, Sand and Teller [41]
iteratively reweight the data term using Gaussian weights. Zis
corresponds to the use of the Welsch loss for in (39). Ze two
mechanisms [41] and [42] are thus equivalent, up to the choice
of the loss function . However, solving (41) is preferable over
(38) since (41) is a smooth optimization problem and it is fully
implicit in —it need no longer be computed. Ze reweighting
method can additionally incorporate weights based on forward
backward flow consistency [43] and flow divergences [41]. Ze
smooth problem (41) moreover admits numerical methods that
are linearly convergent i.e. (1
𝑘
⁄
) whereas problem (38) can
only be solved using numerical methods with a (1
2
⁄
) order
of convergence [45].
Despite the benefits of problem formulation (41), the weight
function associated with the loss (40) would require specifying
a hard threshold, which can cause a discontinuous behavior for
the estimates [36, Sec. 6.4] during the coarse-to-fine estimation
procedure. We use the Welsch loss function in the data term of
(39) since its Gaussian weighting function allows the occlusion
states to vary smoothly from coarser to finer resolutions, and it
does not need hard thresholding. Ze non-convexity introduced
by the use of Welsch loss can be controlled with our graduated
non-convexity algorithm (Section IX).
IX. GN$jTN3OM"()(kT"($!
Ze problem of estimating large motion is nonconvex in two
different senses. First, since a pair of images typically contains
patterns or texture—peaks and troughs—estimating the motion
against these features is generally nonconvex, if the underlying
motion is larger than one pixel. To see this, problem (4) can be
understood as a linearization of the non-convex problem
minimize () =
‖
1
() −
0
‖
2
2
+
2
∗
,
(42)
in which
1
() denotes the sampling of
1
at (+ , ).
Papenberg et al. [46] overcome the data term non-convexity
of (42) with a multiscale (coarse-to-fine) strategy. Zey initially
estimate the motion on coarse or blurred images (and hence the
data term) and the motion is gradually refined by re-estimating
it on gradually finer (less blurred) images. Assuming that there
are scales in total, the motion estimate
𝑗
at the th scale for
≥≥0 can be expressed as the solution of
minimize
𝑗
() =
2
∗

(43)
+
‖
𝑥,𝑗+1
−
𝑗+1
+
𝑡,𝑗+1
‖
2
2
in which
𝑥,𝑗+1
and
𝑡,𝑗+1
denote the - and the -derivatives
of
0
, respectively, evaluated at +
𝑗+1
, , and we assume
the initial flow
0
= .
Ze second type of nonconvexity is due to the particular loss
function used in the data and the regularity terms. Quadratic or
the Huber loss functions could be used so that the estimation is
convex in this second sense, but these convex formulations are
quasi-robust at best and therefore produce flow inappropriately
blurred at flow discontinuities. Zis justifies our use of the non-
convex Welsch loss in the regularity and the data fidelity terms
c0HE3DE3T3K*70LA3GK3H/*6,*2463R4L-.I3LG--3K,+.20G+-E3(+020*LLA:3*3.G+54a0\463LG--3
K,+.20G+3W𝑘 = 0
Z30-3,-463KG/32I43G12070i*20G+:3*+632I43H/*6,*20G+31*/*7424/3
𝑘30-3
0+./4*-4632GJ*/6-3𝐾E3C070L*/LA32G32I432/,+.*2463_,*6/*20.3l9Dm:32I43R4L-.I3LG--3
K,+.20G+3-*2,/*24-32G32I43-*7435*L,4E3CIGJ+3KG/3𝐾 = 2
E
=
0
=
1
=
2
3978 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 29, 2020
of (31) or (39), provided that we solve the resulting problem in
a global manner. Our robust formulation thus has the form
minimize () =
2
‖
(
∗
)
‖
1
1
+
‖
(
𝑥
+
𝑡
)
‖
2
2
(44)
in which denotes the Welsch loss (36).
Our graduated nonconvexity scheme is conceptually similar
to the coarse-to-fine strategy—we initially solve problems (31)
or (39) with a more convex , and refine this initial solution by
solving the same problem using gradually nonconvex (hence
the name graduated non-convexity). Note that Welsch loss (36)
is convex inside the range ± since its second derivative ′′ is
non-negative in this range. We can therefore graduate the non-
convexity of in (31) and (39) by setting sufficiently large at
first, gradually reducing its value to a target filter scale. If there
are graduations, the th graduated motion estimate
𝑘
is the
solution of
minimize
𝑘
() =
2
‖
𝑘
(
∗
)
‖
1
1
(45)
+
‖
𝑘
𝑥
−
𝑗−1
+
𝑡
‖
2
2
in which
𝑘−1
(
,
)
= 2
−2
(
𝐾−𝑘
)
(
, 2
𝐾−𝑘
)
(46)
is a horizontal stretching of the Welsch loss by 2
𝐾−𝑘
. Ze scale
factor 2
−2(𝐾−𝑘)
is needed such that
𝑘
is stretched horizontally
but not vertically. Fig. 6 plots
𝑘
for select . Problem (45) can
be turned into a sequence of reweighted least-squares problems
similar to (43).
To use our graduated nonconvexity scheme together with the
coarse-to-fine estimation strategy, we simply perform one after
the other. Ze overall scheme thus becomes:
procedure graduation
initialize
0
=
(47)
for = , . . . ,1,0
𝑗
←argmin
𝑗
()
end
for = 1,2, . . . ,
𝑘
←argmin
𝑘
()
end
end
in which
𝑗
and
𝑘
are given by (43) and (45) respectively, and
our final motion estimate is
𝐾
. Table I shows the effect of our
graduated optimization scheme on the Middlebury optical flow
dataset. We use = 8 and 2 for the Welsch loss in data and the
regularization terms, respectively. Graduation does not seem to
provide accuracy improvement for Hydrangea, but the value of
its objective (44) does indeed decrease across the GNC stages.
X. N&)P#('TN3OM"()(kT"($!
To solve least-squares problems of the general form (43), we
can use the conjugate-gradient least-squares (CGLS) algorithm
[47]. Let us use the second-order graph Laplacian 
−2
, and
consider the non-stereo flow estimation problem directly. If we
denote = 
𝑥
∗
,
𝑦
∗
∗
, our least-squares problems become
minimize
𝑗+1
() =
2
𝑗
†
𝑗
𝑥
2
2
+
2
𝑗
†
𝑗
𝑦
2
2
(48)
+
(
𝑥,𝑗
,
𝑦,𝑗
)−
𝑗
+
𝑡,𝑗
2
2
,
in which
𝑥,𝑗
,
𝑦,𝑗
,
𝑗
,
𝑗
and
𝑡,𝑗
incorporate the Gaussian
reweighting based on the solution
𝑗
.
In exact arithmetic, problem (48) can be solved using any of
the variants of conjugate gradients. However, in finite precision
arithmetic, the numerical accuracy of the obtained solution can
drop significantly when the usual conjugate gradient method is
used to solve least-squares problems [48, Sec. 7.1]. Conjugate-
gradient least-squares improves the accuracy of the solution by
solving the original, over-constrained system of equations, and
this reduces the impact of round-off errors. Conjugate-gradient
least-squares therefore requires us to express problem (48) as
minimize
𝑗+1
() =
‖
−
‖
2
2
,
(49)
3 33
#,884/RI*L43
UA6/*+H4*3
%0742/G6G+3
h4 + , -3
&/8*+;3
&/8*+93
d/G54;3
d/G5493
T54/* H43
3 J]G3d!'3
3 C2*H43=3
3 C2*H43;3
3 C2*H4393
;EBB3W<E<BVZ3
;E>S3W<E<B<Z3
;EV=3W<E<>nZ3
2.263W0.072Z3
;E;D3W<E=B>Z3
2.073W0.180Z3
;E==3W<E=n;Z3
;E=B3W<E=n>Z3
=E>n3W<E<B;Z3
=ED<3W<E<n9Z3
=ESn3W<E<n;Z3
1.563W0.081Z3
9E9=3W<E;99Z3
9E=V3W<E;;nZ3
9E<V3W<E;;<Z3
3.043W0.219Z3
9E><3W<E9<>Z3
;E>n3W<E;S>Z3
;E9B3W<E;9=Z3
2.203W0.226Z3
9EV;3W<E><;Z3
=EBB3W<E9=BZ3
=EnV3W<E;nBZ3
1.793W0.280Z3
=EBn3W<E=VSZ3
=ES93W<E=<BZ3
=EV;3W<E=<;Z3
1.383W0.099Z3
SE<S3W<EVn9Z3
VEV<3W<EV;<Z3
VE=n3W<E9BDZ3
4.173W0.394Z3
9E<D3W<E;n;Z3
;ES93W<E;==Z3
;E9>3W<E=B>Z3
2.323W0.195Z3
"TjNP3(3
TVERAGE ANGULAR AND WENDQPOINTZ ERROR OF THE GRADUATED NONQCONVEXITY STAGES ON THE )IDDLEBURY TRAINING SET
3 33
#,884/RI*L43
UA6/*+H4*3
%0742/G6G+3
h4 + , -3
&/8*+;3
&/8*+93
d/G54;3
d/G5493
T54/* H43
3 M/G1G-463
3 )61cLGJ
4
3 P10.cLGJ
5
3 'L*--0.!N3
3 !N"h3
3 UC3
2.263W0.072Z3
;ES;3W<E<n;Z3
9E>;3W<E==<Z3
;E9S3W<E<>9Z3
;ESS3W<E<n<Z3
9En<3W<E==nZ3
;E=B3W<E=n>Z3
=EBB3W<E=DVZ3
=EBn3W<E=DVZ3
1.833W<E=S=Z3
=EB<3W0.150Z3
;E;;3W<E=B<Z3
1.563W0.081Z3
9E<;3W<E=S;Z3
=ED;3W<E<n=Z3
;ES>3W<E=9=Z3
9E;=3W<E=D<Z3
VED>3W<E;;VZ3
3.043W0.219Z3
9E;S3W<E;;=Z3
VE;93W<E9<=Z3
9E9V3W<E;9nZ3
VE;>3W<E;B<Z3
SES93W<E99>Z3
;E;<3W<E;;DZ3
1.893W0.181Z3
;E;B3W<E;VBZ3
;E<D3W<E;=BZ3
;EB;3W<E9S<Z3
VE<D3W<EVSBZ3
1.793W0.280Z3
VE>=3W<EVDBZ3
VE=V3W<En=DZ3
;ESn3W<E9>>Z3
9E9D3W<EV9<Z3
>ES;3W<EnSDZ3
1.383W0.099Z3
=En<3W<E=;nZ3
;E9S3W<E=SSZ3
=EVn3W<E=<9Z3
=EBV3W<E=V<Z3
;EnS3W<E;<VZ3
4.17 W0.394Z
VEnn3W<EVD9Z3
DE=B3W<EDS=Z3
SE<;3W<EVD>Z3
DE<S3W<ESV<Z3
DEn=3W<EDB<Z3
2.323W0.195Z3
9E<=3W<E;9;Z3
9E9=3W<E9=DZ3
;EDS3W<E;;<Z3
9E;n3W<E;D>Z3
VEDn3W<E9nSZ3
"TjNP3((3
TVERAGE ANGULAR AND WENDQPOINTZ ERROR OF THE DIFFERENT METHODS ON THE )IDDLEBURY TRAINING SET
3 33
)61cLGJ3
!N"h
6
3
'L*--0.!N3
UC3
P10.cLGJ3
M/G1G-463
3 #,++0+H320743
V<>-3
=S-3
9;9-3
9D-3
=S-3
;9D-3
"TjNP3(((3
PXECUTION TIMES OF THE DIFFERENT METHODS ON3
THE
640 × 480
o&RBAN9p SEQUENCE
4
%4K*,L23)066L48,/A3-4220+H-3,-46E3-7GG2IX3D:3L*/H47G20G+X3<:3G..L,-0G+X3=E3
5
%4K*,L2346H437G64L3*+6364417*2.I0+H3,-463J02I32I43Q7066L48,/A3g*HE
6
dM&32070+H374*-,/463*+63/41G/24638A32I43*,2IG/-3GK3l;mE
3
YOUNG et al.: GRAPH LAPLACIAN REGULARIZATION FOR ROBUST OPTICAL FLOW ESTIMATION 3979
in which
=
𝑥,𝑗
𝑦,𝑗

𝑗
†
𝑗

𝑗
†
𝑗
, =
(
𝑥
,
𝑦
)
𝑗
𝑗
−
𝑡,𝑗
,
(50)
noting that is not explicitly constructed. Only the evaluation
of  and
∗
is required by the conjugate-gradient
least-squares method.
Ze mapping  can be evaluated efficiently during the
optimization since
𝑦,𝑗
and
𝑥,𝑗
are diagonal operators while
the matrix
𝑗
†
𝑗
= −
𝑗
†
𝑗
, where
𝑗
†
is a diagonal matrix
and
𝑗
is a seven-dimensional Gaussian filter. We evaluate the
mapping
𝑗
efficiently using the Gaussian D-tree [40]
filter implementation. Ze mapping
∗
can be evaluated
in a similar way to .
XI. EqMP#()P!"TN3RPC&N"C
For validation of our method, we use the Middlebury optical
flow dataset [49] comprising both natural and synthetic picture
sequences, as well as the newer MPI-Sintel [50] dataset, which
exhibits larger motion. We also compare our method with other
approaches based on variational models. We perform low-level
tasks such as the resampling of the images and flow fields with
built-in MATLAB functions.
A. Middlebury Dataset
Using the Middlebury optical flow dataset, we first study the
effect of filter scales on the accuracy of estimation. Natural and
synthetic image sequences have different characteristics, so we
find the optimal filter parameters using separate grid search for
each sequence type. Ze optimal filter scales are
𝑋,𝑌
= 3 and
𝑍
= 8 for natural sequences, while
𝑋,𝑌
= 6 and
𝑍
= 3 for
synthetic ones. Fig. 7 plots the influence of the two filter scales
on the flow error. For each plot, the other filter scale parameter
is fixed at its optimal value.
Table II lists the average end-point and the angular errors of
the proposed and other well-performing methods on the classic
Middlebury dataset. ClassicNL [1] and MdpFlow [51] methods
are based on continuous optimization similarly to ours whereas
more recent EpicFlow [8] initializes motion estimation with an
interpolation of discrete matches, then refines the estimate with
continuous optimization. Ze non-local total-variation (NLTV)
optical flow estimation method [2] bears some similarities with
our method, but its regularizer is essentially unnormalized, and
uses the Huber loss. Our Laplacian is a fourth-order differential
operator, so it does not penalize the affine trends in the solution
(Section VII). Ze usefulness of preserving affine trends can be
noticed from the results for the Dimetrodon sequence, in which
the motion is piecewise smooth but not translational. Ze affine
interpolator of EpicFlow also preserves affine trends, while the
continuous optimization methods ClassicNL, MdpFlow, NLTV
and Horn-Schunck do not, which incurs higher errors. Table III
provides the execution times of the different method on a CPU
for the 640 ×480 Urban3 sequence. A parallel implementation
of conjugate-gradient least-squares on a GPU can speed up our
method, considering its high potential for parallelization.
To investigate the source of these large improvements in our
estimated flow fields, we include in Fig. 8, a visual comparison
of the motion estimates produced by the different methods and
their corresponding error maps. We observe that improvements
in the motion are mostly near motion discontinuities, as seen in
highly-discontinuous Urban3 and Grove3 sequences. Nonlocal
convex methods like ClassicNL, MdpFlow and NLTV produce
substantially higher errors than our method, suggesting that the
role of our Gaussian edge reweighting (Welsch loss) is equally
important as nonlocal graph regularization itself.
B. MPI-Sintel Dataset
Since large motion estimation methods often use variational
optimization (optical flow estimation) as a post-processing step
to refine the flow [6], [8], we also compare the flow refinement
c0HE3>E3F434e4.23GK3-1*20*L3WL4K2Z3*+63.GLG/3W/0HI2Z3\L24/3-.*L4-3G+32I43*54/*H43*+H,L*/34//G/3*./G--32I43)066L48,/A3-4_,4+.4-E3cG/34*.I31LG2:32I43G2I4/3\L24/3-.*L430-3
\a463*2302-3G1207*L35*L,4E3F43G1207*L3\L24/3-.*L4-3*/43W𝜎
𝑋,𝑌
= 3
:3
𝜎
𝑍
= 8
Z3KG/3+*2,/*L3-4_,4+.4-:3*+63W
𝜎
𝑋,𝑌
= 6
:3
𝜎
𝑍
= 3Z3KG/3-A+2I420.3-4_,4+.4-E3TLL3.,/54-3*/43
\23_,*6/*20.*LLAE
synthetic
natural
synthetic
natural
3980 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 29, 2020
performance of different variational methods on the MPI-Sintel
dataset. We initialize the motion of all methods with the output
from the interpolator of [8]. Using grid search, we find that the
optimal filter parameters
𝑋,𝑌
= 6,
𝑍
= 16 and = 8. Since
the magnitude of the temporal and the flow differences depend
upon the scale of the flow (which can vary considerably across
Sintel sequences), we set the scale parameters of for the data
fidelity and the regularity terms, respectively, as
𝐷
= 0.5 + 0.25,
𝑈
= 0.5 + 0.25,
(51)
in which denotes the standard deviation of the magnitudes of
the initial flow vectors. We use 3 warping stages (1 scale), with
24 conjugate-gradient iterations for each warping stage. We do
not graduate our optimization, but we still allow the ClassicNL
c0HE3nE3h0-,*L0i*20G+3GK32I43/4.G+-2/,.2463gGJ3\4L63W\/-23/GJZ3*+632I434//G/37*13W-4.G+63/GJZ3KG/35*/0G,-307*H43-4_,4+.4-E3j/0HI24/3W/4-1E36*/f4/Z310a4L-3.G//4-1G+63
2G3L*/H4/3W/4-1E3-7*LL4/Z3gGJ354.2G/37*H+02,64-E3F434//G/37*1-3*/43.L011463*23*37*H+02,643GK31
310a4LE
MdpFlow EpicFlow ClassicNL Proposed
Dimetrodon (natural, 583 ×388)
Urban3 (synthetic, 640 ×480)
Grove3 (synthetic, 640 ×480)
YOUNG et al.: GRAPH LAPLACIAN REGULARIZATION FOR ROBUST OPTICAL FLOW ESTIMATION 3981
method to use its own graduated non-convexity scheme.
Table IV provides the average flow accuracy of the different
methods on 14 longer Sintel sequences, each of which contains
50 frames. FullFlow [6] and EpicFlow use the same variational
optimizer, and their errors are shown together. Our method has
lower errors than ClassicNL and the variational optimizer used
by FullFlow and EpicFlow. Note that the Horn-Schunck errors
are occasionally higher than those of the initial flow for certain
sequences.
XII. C$!'N&C($!
In this paper, we analyzed and also proposed an optical flow
estimation approach using dense graph Laplacians. We showed
the advantages of regularizing by dense graph Laplacians from
noise-filtering and edge-preserving perspectives. In the case of
Gaussian edge weighting, we are able to optimize the resulting
objective using fast Gaussian filtering techniques.
Reinterpreting the lagged diffusivity method as re-weighting
of graph edges, we proposed to use the nonconvex Welsch loss
to robustify optical flow estimation against flow discontinuities
and occlusions. Since the weight function of the Welsch loss is
also a Gaussian, our Welsch loss formulation can be optimized
efficiently using fast Gaussian filtering techniques. Our Welsch
formulation produces good flow estimates compared to similar
formulations based on convex losses. A potential disadvantage
of all variational methods is that they involve solving a system
of equations, which is computationally expensive. On the other
hand, variational methods provide a way to accurately estimate
sub-pixel motion. To improve our computational efficiency, we
recommend initializing the flow with a large motion estimation
method, and refining it further using our method.
AMMP!%(q3A
OM"()(k(!d3O&"3"UP3O''N&C($!3VT#(TjNP
To optimize out the occlusion variables from the objective
of (38), first observe that is present only in the distortion and
occlusion terms. Denoting () = 2
‖
log(abs() + 1)
‖
1
1
, let us
define
(
)
= 
𝛽ℎ
(
)
= argmin 
(
)
+
‖
−
‖
2
2
,
(52)
in which the proximal point =
𝑥
+
𝑡
in our context. For
simplified presentation, we assume = 1 since this guarantees
the convexity of the objective of (52) as well. We can write the
first-order optimality conditions of (52) as:
∈(− ) +
‖
log(abs() + 1)
‖
1
1
,
(53)
in which denotes the sub-differential operator. Since problem
(53) is separable across the components of , we now consider
each element of separately.
Examining the first case in which the optimizer
opt
≠0, we
have the relationship
0 =
opt
−+ sgn(
opt
)
1
|
opt
| + 1
,
(54)
implying > +1 if
opt
> 0, and < −1 if
opt
< 0. One can
use the quadratic formula to obtain
opt
=
1
2
−1 +
2
+ 2−3 if > +1,
+ 1 −
√
2
−2−3 if < −1,
(55)
choosing the root away from 0 in both cases. (Ze root close to
0 has the negated sign, which contradicts our assumption about
the sign of
opt
.)
For the remaining case in which
opt
= 0, the sub-derivative
|
log(abs() + 1)
|
=
[
−1, +1
]
, and therefore
opt
= 0 ∈−+
[
−1, +1
]
,
(56)
such that
opt
= 0 if
|
|
≤1. Finally, substituting (55) and (56)
back into the objective of (52) yields (40).
Ze same approach can be used to resolve the proximal map
(52) for all ≤1. When > 1, the objective of (52) becomes
nonconvex with two minima, and the values of the two minima
must be compared to resolve (52). In either case, problem (52)
needs to be solved only once to derive the equivalent loss .
RPcP#P!'PC
3 l=m3 %E3 C,+:3 CE3 #G2I:3 *+63 )E3 @E3 jL*.f:3 oT3 _,*+202*20543 *+*LA-0-3 GK3 .,//4+23
1/*.20.4-30+3G120.*L3gGJ34-207*20G+3*+632I431/0+.01L4-384I0+632I47:p33(/4,
74,"#8'./4,9*&4:35GLE3=<D:3+GE3;:311E3==Sr=9>:3@*+E3;<=VE3
3 l;m3 )E3R4/L84/H4/:3"E3MG.f:3*+63UE3j0-.IGK:3o)G20G+34-207*20G+3J02I3+G+Q
LG.*L3 2G2*L3 5*/0*20G+3 /4H,L*/0i*20G+:p3 0+3 :;<;, 3===, "#8'./%$, 2#>*%/?,
"#(@%$%(>%, #(, "#8'./%$, 9*&*#(, -(), A-//%$(, B%>#+(*/*#(:3 ;<=<:3 11E3
;VDVr;V>=E3
3 l9m3 "E3 j/Ga3 *+63 @E3 )*L0f:3 oN*/H43 60-1L*.474+23 G120.*L3 gGJX3 %4-./012G/3
7*2.I0+H30+3 5*/0*20G+*L3 7G20G+34-207*20G+:p3 3===, C$-(&4, A-//%$(,D(-E4,
F->04,3(/%EE4:35GLE399:3+GE39:311E3S<<rS=9:3)*/E3;<==E3
3 lVm3 'E3N0,:3@E3O,4+:3*+63TE3"G//*L8*:3oC0K23gGJX3%4+-43.G//4-1G+64+.43*./G--3
-.4+4-3*+6302-3*11L0.*20G+-:p33===,C$-(&4,A-//%$(,D(-E4,F->04,3(/%EE4:35GLE3
99:3+GE3S:311E3B>nrBBV:3)*A3;<==E3
3 lSm3 )E3)4+i4:3'E3U401f4:3*+63TE3d40H4/:3o%0-./4243G12070i*20G+3KG/3G120.*L3
gGJ:p30+3G%$8-(,"#(@%$%(>%,#(,A-//%$(,B%>#+(*/*#(:3;<=S:311E3=Dr;nE3
3 lDm3 bE3 'I4+3 *+63 hE3 sGL2,+:3 oc,LL3 gGJX3 $120.*L3 gGJ3 4-207*20G+3 8A3 HLG8*L3
G12070i*20G+3G54/3/4H,L*/3H/06-:p30+3:;<H,3===,"#(@%$%(>%,#(,"#8'./%$,
9*&*#(,-(),A-//%$(,B%>#+(*/*#(,!"9AB6:3N*-3h4H*-:3!h:3&CT:3;<=D:311E3
V><DrV>=VE3
3 l>m3 jE3sE3ME3UG/+3*+63jE3dE3C.I,+.f:3o%424/70+0+H3G120.*L3gGJ:p3D$/*@4,3(/%EE4:3
5GLE3=>:3+GE3=:311E3=nSr;<9:3T,HE3=Bn=E3
3 lnm3 @E3 #45*,6:3 ME3 R40+i*41K4L:3 kE3 U*/.I*G,0:3 *+63 'E3 C.I706:3 oP10.cLGJX3
P6H4Q1/4-4/50+H3 0+24/1GL*20G+3 GK3 .G//4-1G+64+.4-3 KG/3 G120.*L3 gGJ:p3 0+3
:;<I, 3===, "#(@%$%(>%, #(, "#8'./%$, 9*&*#(, -(), A-//%$(, B%>#+(*/*#(,
!"9AB6:3;<=S:311E3==DVr==>;E3
C4_,4+.43
(+020*L3cLGJ
UC
'L*--0.!N
7
c,LL]P10.
8
M/G1G-46
3
*LL4A=3
*LL4A;3
*78,-I>3
8*78GG=3
8*78GG;3
8*+6*H4=3
8*+6*H4;3
.*54V3
7*/f42;3
7G,+2*0+=
-I*7*+;3
-I*7*+93
-L4410+H=3
2471L4;3
=VE93W<En<Z3
>E=3W<E>VZ3
=<E93W<E>VZ3
=VED3W<EnBZ3
=BEV3W=EB>Z3
==E93W=E<<Z3
==E;3W<ED;Z3
nEB3W9EBVZ3
=;E<3W<En;Z3
;SE=3W=E=<Z3
=<ES3W<ES=Z3
BEn3W<ESBZ3
>E93W<EVBZ3
>En3W;ES=Z3
SE<B3W<EVVZ3
9E<>3W<E9nZ3
=9EV3W;EVVZ3
SE9<3W<EV>Z3
=<En3W;E9;Z3
nEVn3W<EB>Z3
DEn<3W<ES;Z3
=9E>3WSEn;Z3
=;E<3W=E=SZ3
9E=S3W<EV>Z3
VEV>3W<E;VZ3
VE<S3W<E;nZ3
=EBB3W<E=9Z3
BES93WVESVZ3
;E>D3W<E;9Z3
=EB=3W<E;DZ3
DEV;3W<ES>Z3
3.523W0.33Z3
>EV>3W1.54Z3
SE=;3W0.58Z3
VED=3W<E;BZ3
7.003W3.50Z3
>E>V3W0.62Z3
;EV93W<EV=Z3
9E993W0.18Z3
;E>>3W<E=nZ3
=ED;3W<E==Z3
4.263W1.99Z3
9E>=3W<E;nZ3
=EnV3W0.24Z3
VE>V3W0.54Z3
VE<B3W<E9BZ3
>EBD3W=ESDZ3
SESS3W<ED=Z3
VE;S3W<E9<Z3
>E9=3W9ES>Z3
>E9B3W<ED9Z3
;E=D3W<E9nZ3
9ED=3W<E;=Z3
;EV93W<E=>Z3
=E;=3W0.08Z3
VE9=3W;E<=Z
2.403W0.22Z3
1.743W<E;DZ3
4.273W0.54Z3
9ESS3W<E9SZ3
7.253W=ESSZ3
5.053W0.58Z3
3.963W0.29Z3
>E;V3W9ES>Z3
7.313W<ED9Z3
1.863W0.38Z3
2.953W0.18Z3
2.233W0.16Z3
1.103W0.08Z3
VE993W;E<=Z3
3
*54/*H43
=;E=3W=E=BZ3
>E;n3W=EVVZ3
VE9S3W0.77Z3
VE993W<E>nZ3
3.953W0.77Z3
3
"TjNP3(h3
TVERAGE ANGULAR AND WENDQPOINTZ ERROR OF THE FLOW AFTER
VA R I AT I O N A L O P T I M IZAT I O N O N T H E )M(QCINTEL DATASET
7
c,LL3W+G+QK*-2Z3-4220+H-3W.L*--0.t+LZ:3=31A/*7063L454L3*+6393J*/10+H3-2*H4-E3
8
%4K*,L23C0+24L3-4220+H-3WQ-0+24LZ:3=31A/*7063L454L3*+6393J*/10+H3-2*H4-E
3
3982 IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 29, 2020
3 lBm3 kE3#4+:3@E3O*+:3jE3!0:3jE3N0,:3qE3O*+H:3*+63UE3kI*:3o&+-,14/50-463%4413
N4*/+0+H3KG/3$120.*L3cLGJ3P-207*20G+:p30+3J*$/?KL*$&/,DDD3,"#(@%$%(>%,
#(,D$/*M>*-E,3(/%EE*+%(>%:3;<=>E3
3l=<m3 CE3 )40-24/:3 @E3 U,/:3 *+63 CE3 #G2I:3 o&+cLGJX3 &+-,14/50-463 N4*/+0+H3 GK3
$120.*L3cLGJ3R02I3*3j060/4.20G+*L3'4+-,-3NG--:p30+3J*$/?K2%>#(),DDD3,
"#(@%$%(>%,#(,D$/*M>*-E,3(/%EE*+%(>%:3;<=nE3
3l==m3 dE3d0L8G*3*+63CE3$-I4/:3o!G+LG.*L3G14/*2G/-3J02I3*11L0.*20G+-32G307*H43
1/G.4--0+H:p3F.E/*&>-E%,F#)%E4,2*8.E4:35GLE3>:3+GE39:311E3=<<Sr=<;n:3!G5E3
;<<nE3
3l=;m3 ME 3s /uI 4+ 8v IL 3* +6 3hE3s GL 2, +:3oP w. 04 +2 3+ G+ LG .* L3/ 4H ,L */0 i* 20G +3KG /3G 12 0. *L 3
gGJ:p30+3"#8'./%$,9*&*#(,N,=""9,:;<::3;<=;:311E39SDr9DBE3
3l=9m3 #E3#*+K2L:3sE3j/4604-:3*+63"E3MG.f:3o!G+QLG.*L32G2*L3H4+4/*L0i4635*/0*20G+3
KG/3G120.*L3gGJ34-207*20G+:p30+3"#8'./%$,9*&*#(,N,=""9,:;<O:3;<=V:311E3
V9BrVSVE3
3l=Vm3 %E3(E3CI,7*+:3CE3sE3!*/*+H:3ME3c/G--*/6:3TE3$/24H*:3*+63ME3h*+64/HI4A+-2:3
oF43 474/H0+H3 \4L63 GK3 -0H+*L3 1/G.4--0+H3 G+3 H/*1I-X3 Pa24+60+H3 I0HIQ
6074+-0G+*L36*2*3*+*LA-0-32G3+42JG/f-3*+63G2I4/30//4H,L*/36G7*0+-:p33===,
2*+(-E,A$#>%&&4,F-+4:35GLE39<:3+GE39:311E3n9rBn:3)*A3;<=9E3
3l=Sm3 @E3 M*+H3 *+63 dE3 'I4,+H:3 od/*1I3 N*1L*.0*+3 /4H,L*/0i*20G+3 KG/3 07*H43
64+G0-0+HX3 T+*LA-0-3 0+3 2I43 .G+20+,G,-3 6G7*0+:p3 3===, C$-(&4, 38-+%,
A$#>%&&4:35GLE3;D:3+GE3V:311E3=>><r=>nS:3T1/E3;<=>E3
3l=Dm3 ME 3 'E3 U *+ -4+ :3 P*&>$%/%,3(Q%$&%,A$#RE%8&1,3(&*+0/,-(),DE+#$*/08&:35GLE3>E3
C(T):3;<=<E3
3l=>m3 UE3k0774/:3TE3j/,I+:3*+63@E3R40.f4/2:3o$120.3gGJ30+3I*/7G+A:p33(/4,74,
"#8'./4,9*&4:35GLE3B9:3+GE39:311E39Dnr9nn:3@,LE3;<==E3
3l=nm3 )E3@E3jL*.f3*+63ME3T+*+6*+:3oF43/G8,-234-207*20G+3GK37,L201L437G20G+-X3
1*/*742/0.3 *+63 104.4J0-4Q-7GG2I3 gGJ3 \4L6-:p3 "#8'./4, 9*&4, 38-+%,
S()%$&/4:35GLE3D9:3+GE3=:311E3>Sr=<V:3@*+E3=BBDE3
3l=Bm3 cE3 #E3 sE3 'I,+H3 *+63 cE3 'E3 d/*I*7:3 2'%>/$-E, G$-'0, J%#$?E3 T74/0.*+3
)*2I47*20.*L3CG.E:3=BB>E3
3l;<m3 !E3j0HH-:3!E3NE3j0HH-:3*+63jE3!G/7*+:3DE+%R$-*>,G$-'0,J%#$?:35GLE3D>E3
'*78/06H43,+054/-02A31/4--:3=BB9E3
3l;=m3 )E3 h4224/L0:3 @E3 sG5*x450y:3 *+63 hE3 sE3 dGA*L:3 L#.()-/*#(&, #@, 2*+(-E,
A$#>%&&*(+E3'*78/06H43&+054/-02A3M/4--:3;<=VE3
3l;;m3'E3"G7*-03*+63#E3)*+6,.I0:3oj0L*24/*L3\L24/0+H3KG/3H/*A3*+63.GLG/307*H4-:p3
0+32*T/0,3(/%$(-/*#(-E,"#(@%$%(>%,#(,"#8'./%$,9*&*#(:3=BBn:311E3n9BrnVDE3
3l;9m3 TE3j,*64-:3jE3'GLL:3*+63@E3)E3)G/4L:3oT3+G+QLG.*L3*LHG/02I73KG/307*H43
64+G0-0+H:p3 0+3 :;;I, 3===, "#8'./%$, 2#>*%/?, "#(@%$%(>%, #(, "#8'./%$,
9*&*#(,-(), A-//%$(, B%>#+(*/*#(, !"9ABU;I6:3;<<S:35GLE3;:311E3D<rDS35GLE3
;E3
3l;Vm3 @E3PE3%4++0-3@/E3*+63#E3PE3R4L-.I:3o"4.I+0_,4-3KG/3+G+L0+4*/3L4*-23-_,*/4-3
*+63/G8,-23/4H/4--0G+:p3"#88.(4,2/-/4,K,2*8.E4,"#8'./4:35GLE3>:3+GE3V:311E3
9VSr9SB:3@*+E3=B>nE3
3l;Sm3 )E3T8/*7GJ02i3 *+63 (E3C24H,+:3oU*+68GGf3GK3 7*2I47*20.*L3K,+.20G+-:p3
D84,74,A0?&4:35GLE39V:3+GE3;:311E3=>>r=>>:3=BDDE3
3l;Dm3 )E3TE3M0+-fA:33(/$#).>/*#(,/#,L#.$*%$,D(-E?&*&,-(),V-Q%E%/&E3T74/0.*+3
)*2I47*20.*L3CG.E:3;<<nE3
3l;>m3 UE3!*H4L3*+63RE3P+f4L7*++:3oT+30+54-20H*20G+3GK3-7GG2I+4--3.G+-2/*0+2-3
KG/32I434-207*20G+3GK360-1L*.474+2354.2G/3\4L6-3K/G7307*H43-4_,4+.4-:p3
3===,C$-(&4,A-//%$(,D(-E4,F->04,3(/%EE4:35GLE3n:3+GE3S:311E3SDSrSB9:3C41E3
=BnDE3
3l;nm3 )E3 TE3 C+A64/:3 o$+3 2I43 7*2I47*20.*L3 KG,+6*20G+-3 GK3 -7GG2I+4--3
.G+-2/*0+2-3 KG/3 2I43 6424/70+*20G+3 GK3 G120.*L3 gGJ3 *+63 KG/3 -,/K*.43
/4.G+-2/,.20G+:p33===,C$-(&4,A-//%$(,D(-E4,F->04,3(/%EE4:35GLE3=9:3+GE3==:3
11E3==<Sr===V:3!G5E3=BB=E3
3l;Bm3 'E3C.I+z//:3oC4H74+2*20G+3GK350-,*L37G20G+38A370+070i0+H3.G+54a3+G+Q
_,*6/*20.3K,+.20G+*L-:p30+3A$#>%%)*(+&,#@,<:/0,3(/%$(-/*#(-E,"#(@%$%(>%,
#(,A-//%$(,B%>#+(*/*#(:3=BBV:35GLE3=:311E3DD=rDD935GLE=E3
3l9<m3 @E3 R40.f4/23 *+63 'E3 C.I+z//:3 oT3 2I4G/420.*L3 K/*74JG/f3 KG/3 .G+54a3
/4H,L*/0i4/-30+3M%PQ8*-463.G71,2*20G+3GK307*H437G20G+:p33(/4,74,"#8'./4,
9*&4:35GLE3VS:3+GE39:311E3;VSr;DV:3%4.E3;<<=E3
3l9=m3 dE3T,84/2:3#E3%4/0.I4:3*+63ME3sG/+1/G8-2:3o'G71,20+H3G120.*L3gGJ350*3
5*/0*20G+*L324.I+0_,4-:p323DF,74,D''E4,F-/04:35GLE3D<:3+GE3=:311E3=SDr=n;:3
@*+E3=BBBE3
3l9;m3 dE3 C2/*+H:3 "#8'./-/*#(-E, 2>*%(>%, -(), =(+*(%%$*(+E3 R4LL4-L4AQ
'*78/06H43M/4--:3;<<>E3
3l99m3 #E3 )4//0-:3 o%G,8LA3 -2G.I*-20.3 H/*1I3 7*2/0.4-:p3A.RE4, =E%W/$#/%0(*XW#+,
L-W4,2%$4,F-/4:3+GE3n:311E3DVr>=:3=BB>E3
3l9Vm3 RE3 UE3 M/4-- : 3 C E 3 TE3 "4,fGL-fA:3 RE3 "E3 h4224/L0+H:3 *+63 j E 3 ME3 cL*++4/A:3
Y.8%$*>-E, B%>*'%&, Z$), =)*/*#(1, J%, D$/, #@, 2>*%(/*M>, "#8'./*(+E3
'*78/06H43&+054/-02A3M/4--:3;<<>E3
3l9Sm3 "E3 cE3 'I*+3 *+63 @E3 W@*.f04Z3 CI4+:3 38-+%, A$#>%&&*(+, -(), D(-E?&*&1,
9-$ *- /* #(- E[ ,AP =[ ,V- Q% E% /[, -( ), 2/# >0 -& /*> ,F %/0 #) &E3C(T):3;<<SE3
3l9Dm3 RE3 #4 A:3 3(/$#).>/*#(, /#, B#R.&/, -(), \.-&*KB#R.&/, 2/-/*&/*>-E,F%/0#)&E3
j4/L0+3U4064L84/HX3C1/0+H4/Qh4/ L* H:3 =B n9E 3
3l9>m3 kE3kI*+H:3oM*/*7424/34-207*20G+324.I+0_,4-X3*32,2G/0*L3J02I3*11L0.*20G+32G3
.G+0.3\220+H:p338-+%,9*&4,"#8'./4:35GLE3=S:3+GE3=:311E3SBr>D:3@*+E3=BB>E3
3l9nm3TE3jL*f43*+63TE3k0--4/7*+:39*&.-E,B%>#(&/$.>/*#(E3'*78/06H4:3)T:3&CTX3
)("3M/4--:3=Bn>E3
3l9Bm3 sEQ@E3 OGG+3 *+63 (E3 CE3 sJ4G+:3 oT6*120543 -,11G/2QJ40HI23 *11/G*.I3 KG/3
.G//4-1G+64+.43-4*/.I:p33===,C$-(&4,A-//%$(,D(-E4,F->04,3(/%EE4:35GLE3;n:3
+GE3V:311E3DS<rDSD:3T1/E3;<<DE3
3lV<m3 TE3T6*7-:3!E3d4LK*+6:3@E3%GL-G+:3*+63)E3N45GA:3od*,--0*+3s%Q2/44-3KG/3
K*-23I0HIQ6074+-0G+*L3\L24/0+H:p30+3D"F,23GGBDA],:;;^,A-'%$&:3!4J3
OG/f :3! O:3 &CT :3; <<B : 31 1E3;= X=r;=X=;E3
3lV=m3 ME 3C* +6 3* +6 3CE 3"4 LL 4/: 3oM */ 20 .L4 350 64 GX3 LG+ HQ/*+H437G20G+34-207*20G+3,-0+H3
1G0+232/*[4.2G/04-:p33(/4,74,"#8'./4,9*&4:35GLE3n<:3+GE3=:31E3>;:3$.2E3;<<nE3
3lV;m3 TE3TA5*.0:3)E3 #*120-:3 *+63 CE3 CG*22G:3 oC1*/-43 G..L,-0G+3 6424.20G+3 J02I3
G120.*L3gGJ:p33(/4,74,"#8'./4,9*&4:35GLE3B>:3+GE39:311E39;;r99n:3)*A3;<=;E3
3lV9m3 CE3(+.43*+63@E3sG+/*6:3o$..L,-0G+Q*J*/43G120.*L3gGJ34-207*20G+:p33===,
C$-(&4,38-+%,A$#>%&&4:35GLE3=>:3+GE3n:311E3=VV9r=VS=:3T,HE3;<<nE3
3lVVm3 PE3 @E3 '*+6{-:3 )E3 jE3 R*f0+:3 *+63 CE3 ME3 jGA6:3 oP+I*+.0+H3 -1*/-02A3 8A3
/4J40HI2463|=370+070i*20G+:p374,L#.$*%$,D(-E4,D''E4:35GLE3=V:3+GE3S:311E3
n>>rB<S:3%4.E3;<<nE3
3lVSm3 OE 3 ! 4 - 2 4 / G 5: 3 o T3 7 4 2 I G 6 3 K G / 3 - G L 5 0 + H 3 2 I 4 3 . G + 5 4 a 3 1 / G H / * 7 7 0 + H 3 1 / G 8 L 4 7 3
J02I3.G+54/H4+.43/*243$W=]f
2
Z:p3P#WE,DW-),Y-.W,222B:35GLE3;DB:311E3SV9r
SV>:3=Bn9E3
3lVDm3 !E3 M*14+84/H:3TE3 j/,I+:3 "E3 j/Ga:3 CE3 %06*-:3 *+63 @E3 R40.f4/2:3 oU0HILA3
*..,/*243G120.3gGJ3.G71,2*20G+3J02I32I4G/420.*LLA3[,-20\463J*/10+H:p33(/4,
74,"#8'./4,9*&4:35GLE3D>:3+GE3;:311E3=V=r=Sn:3T1/E3;<<DE3
3lV>m3 )E3 #E3 U4-24+4-3 *+63 PE3 C204K4L:3 o)42IG6-3 GK3 .G+[,H*243 H/*604+2-3 KG/3
-GL50+H3L0+4*/3-A-247-:p374,B%&4,Y-/E4,_.$4,2/-()4:35GLE3VB:3+GE3D:311E3V<Br
V9D:3=BS;E3
3lVnm3 'E3'E3M*0H43*+63)E3TE3C*,+64/-:3oNCb#X3T+3*LHG/02I73KG/3-1*/-43L0+4*/3
4_,*20G+-3*+63-1*/-43L4*-23-_,*/4-:p3D"F,C$-(&,F-/0,2#@/`:35GLE3n:3+GE3=:3
11E3V9r>=:3)*/E3=Bn;E3
3lVBm3 CE3j*f4/:3%E3C.I*/-240+:3@E3ME3N4J0-:3CE3#G2I:3)E3@E3jL*.f:3*+63#E3Ci4L0-f0:3
oT36*2*8*-43*+6345*L,*20G+3742IG6GLGHA3KG/3G120.*L3gGJ:p33(/4,74,"#8'./4,
9*&4:35GLE3B;:3+GE3=:311E3=r9=:3)*/E3;<==E3
3lS<m3 %E3@E3j,2L4/:3@E3R,Le:3dE3jE3C2*+L4A:3*+63)E3@E3jL*.f:3oT3+*2,/*L0-20.3G14+3
-G,/.437G5043KG/3G120.*L3gGJ345*L,*20G+:p30+3"#8'./%$,9*&*#(,N,=""9,
:;<::3;<=;:311E3D==rD;SE3
3lS=m3 NE3q,:3@E3@0*:3*+63OE3)*2-,-I02*:3o)G20G+3642*0L31/4-4/50+H3G120.*L3gGJ3
4-207*20G+:p3 0+3:;<;, 3===, "#8'./%$, 2#>*%/?, "#(@%$%(>%,#(, "#8'./%$,
9*&*#(,-(),A-//%$(,B%>#+(*/*#(:3;<=<:311E3=;B9r=9<<E
Sean I. Young3/4.405463jE'G73*+63jEPE364H/44-30+3-GK2J*/434+H0+44/0+H3K/G73
2I43&+054/-02A3GK3T,.fL*+630+3;<<n:3*+632I43)EP+HC.3*+63MIE%E364H/44-3K/G73
2I43&+054/-02A3GK3!4J3CG,2I3R*L4-:30+3;<==3*+63;<=n:3/4-14.2054LAE3',//4+2LA:3
I430-3*31G-26G.2G/*L3/4-4*/.I4/3*23C2*+KG/63&+054/-02A:3C2*+KG/6:3'TE3(+3;<=D:3
I43J*-3*350-020+H3/4-4*/.I4/3J02I3(+24/%0H02*L3'G77,+0.*20G+-:3C*+3%04HG:3'TE3
U0-3 /4-4*/.I3 0+24/4-2-3 */43 L*/H4Q-.*L43 G12070i*20G+3 *+63 0+54/-43 1/G8L47-3 0+3
07*H43 1/G.4--0+HE3 U43 /4.405463 2I43TM#C](TM#3 84-23 1*14/3 *J*/63 *23 %('"T3
;<=n:32GH42I4/3J02I3%*5063"*,87*+E3
Aous T. Naman3/4.405463jEC.E364H/4430+34L4.2/G+0.-3*+6324L4.G77,+0.*20G+3
4+H0+44/0+H3K/G73TLQ!*I/*0+3&+054/-02A:3j*HI6*6:3(/*_:30+3=BBV:32I43)EP+HC.E3
64H/4430+34+H0+44/0+H3K/G732I43&+054/-02A3GK3)*L*A*:3s,*L*3N,71,/:3)*L*A-0*:3
0+3;<<<:3*+632I43MIE%E364H/4430+34L4.2/0.*L34+H0+44/0+H3K/G73F43&+054/-02A3GK3
!4J3CG,2I3R*L4-3W&!CRZ:3CA6+4A:3T,-2/*L0*:30+3;<==E3U430-3.,//4+2LA3*3MG-2Q
6G.2G/*L3#4-4*/.I4/3 J02I3 2I43 C.IGGL3 GK3PL4.2/0.*L3 P+H0+44/0+H3*+63"4L4.G7Q
7,+0.*20G+-:3&!CRE3U0-3/4-4*/.I30+24/4-2-3*/430+307*H43G/35064G3.G71/4--0G+3
*+6364L054/AE3
David Taubman3/4.405463jECE3*+63jEPE364H/44-30+34L4.2/0.*L34+H0+44/0+H3K/G73
2I43&+054/-02A3GK3CA6+4A:30+3=BnD3*+63=Bnn:3/4-14.2054LA:3*+632I43)ECE3*+63MIE%E3
64H/44-3 K/G73 2I43 &+054/-02A3 GK3 '*L0KG/+0*3 *23 j4/f4L4A:3 0+3 =BB;3 *+63 =BBV:3
/4-14.2054LAE3 c/G73 =BBV3 2G3 =BBn:3 I43 J*-3 J02I3 U4JL422QM*.f*/6}-3 #4-4*/.I3
N*8G/*2G/04-:3M*LG3TL2G:3'TE3U43[G0+463&!CR30+3=BBn:3JI4/43I430-3.,//4+2LA3*3
M/GK4--G/3J02I32I43C.IGGL3GK3PL4.2/0.*L3P+H0+44/0+H3*+63"4L4.G77,+0.*20G+-E3
U43 I*-3 *,2IG/463 2I43 8GGf3 @MPd;<<<X3 (7*H43 'G71/4--0G+3 c,+6*74+2*L-:3
C2*+6*/6-3*+63M/*.20.4:3J02I3)E3)*/.4LL0+E3U0-3/4-4*/.I30+24/4-2-30+.L,643I0HILA3
-.*L*8L43 07*H43 *+63 5064G3 .G71/4--0G+:3 7G20G+3 4-207*20G+3 *+63 7G64L0+H:3
0+54/-431/G8L47-30+307*H0+H:314/.412,*L37G64L0+H:3*+637,L207460*360-2/08,20G+3
-A-247-E3U43/4.4054632I43&+054/-02A3)46*L3K/G732I43&+054/-02A3GK3CA6+4AE3U43
I*-3/4.4054632JG384-231*14/3*J*/6-3K/G732I43(PPP3'0/.,02-3*+63CA-247-3CG.042A3
KG/3 2I43 =BBD3 1*14/3 4+202L463T3'G77G+3 c/*74JG/f3 KG/3 #*243 *+63 %0-2G/20G+Q
j*-463C.*L0+H3GK3U0HILA3C.*L*8L43'G71/4--463h064G:3*+63K/G732I43(PPP3C0H+*L3
M/G.4--0+H3 CG.042A3 KG/3 2I43 ;<<<3 1*14/3 4+202L463 U0HI3 M4/KG/7*+.43 C.*L*8L43
(7*H43'G71/4--0G+3J02I3Pj'$"E