IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2017

1

Continuous 3D Label Stereo Matching using
Local Expansion Moves
Tatsunori Taniai, Member, IEEE, Yasuyuki Matsushita, Senior Member, IEEE, Yoichi Sato, Member, IEEE,and Takeshi Naemura, Member, IEEE
Abstract—We present an accurate stereo matching method using local expansion moves based on graph cuts. This new
move-making scheme is used to efficiently infer per-pixel 3D plane labels on a pairwise Markov random field (MRF) that effectively
combines recently proposed slanted patch matching and curvature regularization terms. The local expansion moves are presented as
many α-expansions defined for small grid regions. The local expansion moves extend traditional expansion moves by two ways:
localization and spatial propagation. By localization, we use different candidate α-labels according to the locations of local
α-expansions. By spatial propagation, we design our local α-expansions to propagate currently assigned labels for nearby regions.
With this localization and spatial propagation, our method can efficiently infer MRF models with a continuous label space using
randomized search. Our method has several advantages over previous approaches that are based on fusion moves or belief
propagation; it produces submodular moves deriving a subproblem optimality; it helps find good, smooth, piecewise linear disparity
maps; it is suitable for parallelization; it can use cost-volume filtering techniques for accelerating the matching cost computations. Even
using a simple pairwise MRF, our method is shown to have best performance in the Middlebury stereo benchmark V2 and V3.
Index Terms—Stereo Vision, 3D Reconstruction, Graph Cuts, Markov Random Fields, Discrete-Continuous Optimization.

F

1

I NTRODUCTIONSTEREO vision often struggles with a bias toward reconstructing fronto-parallel surfaces, which can stemfrom matching cost, smoothness regularization, or eveninference [5], [45]. Segment-based stereo [4] that represents
disparity maps by disparity plane segments, could easethe bias issue but recovered surfaces are constrained to bepiecewise planar. Recently, two breakthroughs have been
independently made that overcome the fronto-parallel bias
while surpassing the limitation of segment-based stereo.
One breakthrough is matching cost using slanted patch
matching [5]. In this approach, the disparity dp of each pixel
p is over-parameterized by a local disparity plane

dp = ap u + bp v + cp

(1)

defined on the image domain (u, v), and the triplet(ap , bp , cp ) is estimated for each pixel p instead of directly
estimating dp . The matching window is then slanted according to Eq. (1), which produces linearly-varying disparities
within the window and thereby measures patch dissimilarity accurately even with large matching windows.
Another key invention is curvature regularization [27],[33] by tangent plane labels. Unlike second-order smoothness [45] that forms a higher-order term, this curvature
regularization is nicely represented by pairwise terms likeconventional (fronto-parallel) linear [19] and truncated linear models, and can handle smooth surfaces beyond planes.
Given slanted patch matching [5] and tangent-based
curvature regularization [33], the use of 3D plane labels
allows us to establish a new stereo model using a pairwise
Markov random field (MRF) [12] that is free from the••T. Taniai, T. Naemura, and Y. Sato are with the University of Tokyo, Japan.
Y. Matsushita is with Osaka University, Japan.

Manuscript received Apr. 22, 2015; revised Jul. 4 and Oct. 12, 2017.

Random init.

1 iteration

2 iterations

Post-processing

Fig. 1. Evolution of our stereo matching estimates. From top to bottom, we show disparity maps, normal maps of disparity planes, anderror maps with 0.5 pixel threshold where ground truth is given. In ourframework, we start with random disparities that are represented by
per-pixel 3D planes, i.e., disparities (top) and normals (middle). We
then iteratively apply our local expansion moves using GC (middles)
to update and propagate local disparity planes. Finally, the resulting
disparity map is further refined at a post-processing stage using leftright consistency check and weighted median filtering (rightmost).

fronto-parallel bias. However, while stereo with standard
1D discrete disparity labels [8], [23], [42] can be directly
solved by discrete optimizers such as graph cuts (GC) [7],[24] and belief propagation (BP) [11], [47], such approaches
cannot be directly used for continuous 3D labels due to thehuge or infinite label space (a, b, c) ∈ R3 .
To efficiently infer 3D labels, recent successful methods [3], [5], [14], [32] use PatchMatch [1], [2], an inference
algorithm using spatial label propagation. In PatchMatch,each pixel is updated in raster-scan order and its refinedIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2017label is propagated to next pixels as their candidate labels.
Further in [3], this sequential algorithm is combined withBP yielding an efficient optimizer PMBP for pairwise MRFs.
In terms of MRF optimization, however, BP is considereda sequential optimizer, which improves each variable individually keeping others conditioned at the current state.
In contrast, GC improves all variables simultaneously byaccounting for interactions across variables, and this global
property helps optimization avoid bad local minimums [38],[45]. In order to take advantage of this and efficiently infer
3D planes by GC, it is important to use spatial propagation. Nevertheless, incorporating such spatial propagation
into GC-based optimization is not straightforward, because
inference using GC rather processes all nodes simultaneously,
not one-by-one sequentially like PatchMatch and BP.
In this paper, we introduce a new move making scheme,
local expansion moves, that enables spatial propagation in GC
optimization. The local expansion moves are presented as
many α-expansions [8] defined for a small extent of regions
at different locations. Each of this small or local α-expansion
tries to improve the current labels in its local region in anenergy minimization manner using GC. Here, those current
labels are allowed to move to a candidate label α, which is
given uniquely to each local α-expansion in order to achieve
efficient label searching. At the same time, this procedure isdesigned to propagate a current label in a local region for
nearby pixels. For natural scenes that often exhibit locallyplanar structures, the joint use of local expansion moves
and GC has a useful property. It allows multiple pixels ina local region to be assigned the same disparity plane by asingle min-cut in order to find a smooth solution. Being able
to simultaneously update multiple variables also helps to
avoid being trapped at a bad local minimum.
For continuous MRF inference, fusion-based approaches
using GC [25] are often employed with some heuristics togenerate disparity map proposals [33], [34], [45]. Our work
bridges between apparently different inference approaches
of PatchMatch and GC, and can take advantage of thoseheuristics used in both PatchMatch and fusion based approaches leading to higher efficiency and accuracy.
The advantages of our method are as follows. 1) Our
local expansion move method produces submodular moves
that guarantee the optimal labeling at each min-cut (subproblem optimal), which in contrast is not guaranteed ingeneral fusion moves [25]. 2) This optimality and spatial
propagation allow randomized search, rather than employ
external methods to generate plausible initial proposals as
done in fusion approaches [25], [45], which may limit thepossible solutions. 3) Our method achieves greater accuracy
than BP [3] thanks to the good properties of GC and local
expansion moves. 4) Unlike other PatchMatch based methods [3], [5], [14], our method can incorporate the fast costfiltering technique of [32]. In this manner, we can efficiently
reduce the computation complexity of unary terms from
O(|W |) to approximately O(1), removing dependency fromsupport window size |W |. 5) Unlike PMBP [3], our method
is well suited for parallelization using both CPU and GPU.1
1. Although BP is usually parallelizable, PMBP differs from BP’s
standard settings in that it defines label space dynamically for each pixel
and propagates it; both make parallelization indeed non-trivial.

2

With multiple CPU cores, each of our local α-expansions(i.e., min-cut computations) can be individually performedin parallel. With a GPU implementation, computations ofunary terms can be more efficiently performed in a parallelmanner for further acceleration.
This paper is an extended version of our conference
paper [39]. The extensions are summarized as follows. We
add theoretical verifications on the preferability of our
method for piecewise planar scenes in Sec. 3.1, and also
on the subproblem optimality of our algorithm in Sec. 3.3.
Furthermore, the efficiency of our algorithm is improved by
two ways; In Sec. 3.3 we show the parallelizablity of our
local expansion move algorithm; In Sec. 3.5 we show thatthe fast cost filtering technique of [32] can be used in our
method. The effectiveness of both extensions is thoroughly
evaluated in the experiments, and we show that even a CPU
implementation of the proposed method achieves about
2.1x faster running times than our previous GPU implementation [39], with comparable or even greater accuracy.
Finally, we add experimental evaluations in Sec. 4 showingthat our method outperforms a fusion-based approach [33]and all the state-of-the-art methods registered in the latest
Middlebury benchmark V3. Our code is publicly available
online (https://github.com/t-taniai/LocalExpStereo).

2

R ELATED WORK

2.1

MRF stereo methods

MRF stereo methods can be categorized into thee approaches: discrete stereo, segment-based stereo, and continuous stereo.
Discrete stereo [8], [23], [42] formulates stereo matchingas a discrete multi-labeling problem, where each pixel isindividually assigned one of pre-defined discrete disparity
values. For this problem, many powerful discrete optimizers, such as BP [11], [47], TRW [21], and GC [7], [24], can bedirectly used. Successful results are shown using GC withexpansion moves [8], [38]. In expansion moves, the multilabeling problem is reduced to a sequence of binary-labelingproblems, each of which can be exactly solved by GC, if only
pairwise potentials ψ meet the following submodularity of
expansion moves [8], [22]:ψ(α, α) + ψ(β, γ) ≤ ψ(β, α) + ψ(α, γ).

(2)Segment-based stereo [4], [15], [41], [43] assigns a 3D
disparity plane for each of over-segmented image regions.
Although this approach yields continuous-valued disparities, it strictly limits the reconstruction to a piecewise planar
representation and is subject to the quality of initial segmentation. More recent methods alleviate these limitations
by using a complex layered MRF [6], multi-scale segmentation [26], or jointly estimating segmentation [46].
The last group, to which our method belongs, is continuous stereo [3], [5], [14], [32], [33], [45], where each pixel
is assigned a distinct continuous disparity value. Somemethods [33], [45] use fusion moves [25], an operation thatcombines two solutions to make better one (binary fusion)
by solving a non-submodular binary-labeling problem using
QPBO-GC [22]. In this approach, a number of continuousvalued disparity maps (or so-called proposals in the literature [25]) are first generated by other external methods (e.g.,

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2017segment-based stereo [45]), which are then combined as asequence of binary fusions. Our method differs from thisfusion-based approach in that we use spatial propagation
and randomization search during inference, by which we
only require a randomized initial solution instead of thosegenerated by external methods. Also, binary energies produced in our method are always submodular allowing exact
inference via GC (subproblem optimal). More importantly,
we rather provide an efficient inference mechanism that isconceptually orthogonal to proposal generation. Thus, conventional proposal generation schemes used in fusion approaches [33], [34], [45] can be incorporated in our method
for further improvements. A stereo method by Bleyer et
al. [5] proposes accurate photo-consistency measures using
3D disparity planes that are inferred by PatchMatch [1],[2]. Heise et al. [14] incorporate Huber regularization into
[5] using convex optimization. Besse et al. [3] point out aclose relationship between PatchMatch and BP and present
a unified method called PatchMatch BP (PMBP) for pairwise
continuous MRFs. PMBP is probably the closest approach toours in spirit, but we use GC instead of BP for the inference.
Therefore, our method is able to take advantage of better
convergence of GC [38] for achieving greater accuracy. In
addition, our method allows efficient parallel computationof unary matching costs and even min-cut operations.
2.2

3

their subregion filtering technique can be effectively incorporated into our method, and we achieve greater accuracy
than their local stereo method [32].

3

P ROPOSED METHODThis section describes our proposed method. Given two
input images IL and IR , our purpose is to estimate thedisparities of both images.
In Sec. 3.1, we first discuss a geometric property of ourmodel. We then define our energy function in Sec. 3.2, anddescribe the fundamental idea of our optimizing strategy
and its properties in Sec. 3.3. The whole optimization procedure is presented in Sec. 3.4, and we further discuss a fast
implementation in Sec. 3.5.
3.1

Geometric interpretation to slanted patch matching

Slanted patch matching [5] by Eq. (1) implicitly assumes thatthe true disparity maps are approximately piecewise linear.
While it is not discussed in [5], linear disparity is exactly
related to motion by planar surfaces [4], [33]. Formally, if there
exists a plane in the 3D world coordinates (x, y, z) ∈ R3

a0p x + b0p y + c0p z = h0p

parameterized by (a0p , b0p , c0p , h0p ), then motion due to thisgeometric plane is represented by a disparity plane

Cost-volume filtering

Patch-based stereo methods often use cost-volume filtering
for fast implementations. Generally, computing a matching
cost C for a patch requires O(|W |) of computation, where|W | is the size of the patch. However, given a cost-volume
slice ρd (p) that represents pixelwise raw matching costs
kI(p) − I 0 (p − d)k for a certain disparity label d, the patchbased matching costs can be efficiently computedby apP
plying a filtering to the cost map as Cd (p) = q ωpq ρd (q).
Here, the filter kernel ωpq represents the matching window
at p. If we use a constant-time filtering, each matching cost
Cd (p) is efficiently computed in O(1).
The box filtering can achieve O(1) by using integral
image but such simple filtering flattens object boundaries.
For this boundary issue, Yoon and Kweon [48] propose anadaptive support-window technique that uses the joint bilateral filtering [35] for cost filtering. Although this adaptive
window technique successfully deals with the boundary
issue [16], it involves O(|W |) of computation because of thecomplexity of the bilateral filtering. Recently, He et al. [13]propose a constant-time edge-aware filtering named theguided image filtering. This filtering is employed in a costvolume filtering method of [17], [36], achieving both edgeawareness and O(1) of matching-cost computation.
In principle, stereo methods using PatchMatch inference [3], [5], [14] cannot take advantage of the cost filtering
acceleration, since in those methods the candidate disparity
labels are given dynamically to each pixel and we cannot
make a consistent-label cost-volume slice ρd (p). To thisissue, Lu et al. [32] extend PatchMatch [5] to use superpixels as a unit of cost calculations. In their method, called
PatchMatch filter (PMF), fast cost-volume filtering of [17],[36] is applied in small subregions and that approximately
achieves O(1) of complexity. We will show in Sec. 3.5 that(3)d(u, v) =B 0
ap u + b0p v + f c0p ,
0
hp

(4)

where B and f are the baseline and focal length of the stereo
cameras.2 This piecewise planar assumption is reasonable
as many natural scenes approximately exhibit locally planar
surfaces. As we will discuss later, this fact also supports our
method design as our model and inference both softly preferpiecewise linear disparity maps while having the ability to
express arbitrary curved surfaces.
3.2

Formulation

We use the 3D plane label formulation to take advantageof the powerful slanted patch matching [5] and curvature
regularization [33]. Here, each pixel p’s disparity dp is overparameterized by a 3D plane dp = ap u + bp v + cp . Therefore, we seek a mapping fp = f (p) : Ω → L that assigns a
disparity plane fp = (ap , bp , cp ) ∈ L for every pixel p in theleft and right images. Following conventional MRF stereo
methods [8], [23], [33], [42], we estimate f by minimizing
the following energy function based on a pairwise MRF.
X
X
E(f ) =
φp (fp ) + λψpq (fp , fq ).
(5)p∈Ω(p,q)∈N

The first term, called the data term or unary term, measures
the photo-consistency between matching pixels. The disparity plane fp defines a warp from a pixel p in one image toits correspondence in the other image. The second term iscalled the smoothness term or pairwise term, which penalizes
discontinuity of disparities between neighboring pixel pairs(p, q) ∈ N . We define these terms below.
2. Derived easily by using the perspective projection x = uz/f andy = vz/f , and the depth-disparity relation z = Bf /d in rectified stereo.

���������

Data term
To measure photo-consistencies, we use the slanted patch
matching term that has been recently proposed by [5]. The
data term of p in the left image is defined asX
φp (fp ) =ωps ρ(s|fp ).
(6)4���������

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2017

��

��

��

��

�

s∈Wp

�

�

�
�

�

Here, Wp is a square window centered at p. The weight ωps
implements the adaptive support window proposed in [48].
For this we replace the bilateral filter weight used in [5], [48]with more advanced guided image filtering [13] as below.

X 
1
T
−1
ωps =
1
+
(I
−
µ
)
(Σ
+
e)
(I
−
µ
)
p
k
k
s
k
|W 0 |20

Fig. 2. Illustration of the smoothness term proposed in [33]. (a) The
smoothness term penalizes the deviations of neighboring disparity
planes shown as red arrows. (b) When neighboring pixels are assigned
the same disparity plane, it gives no penalty; thus, it enforces second
order smoothness for the disparity maps.

(7)Here, Ip = IL (p)/255 is a normalized color vector, µk andΣk are the mean and co-variance matrix of Ip in a local
regression window Wk0 , and e is an identity matrix with asmall positive coefficient for avoiding over-fitting. We will
discuss how to efficiently compute this filtering in Sec. 3.5.
Given a disparity plane fp = (ap , bp , cp ), the function
ρ(s|fp ) in Eq. (6) measures the pixel dissimilarity between
a support pixel s = (su , sv ) in the window Wp and its
matching point in the right image

Notice that ψ̄pq (fp , fq ) = 2|cp − cq | when a = b = 0
is forced; therefore, the smoothness term ψpq (fp , fq ) naturally extends the traditional truncated linear model [8],although the latter has a fronto-parallel bias and should be
avoided [33], [45]. Also, as shown in Fig. 2 (b), this term
becomes zero when fp = fq . This enforces piecewise linear
disparity maps and so piecewise planar object surfaces.
Furthermore, this term satisfies the following property fortaking advantage of GC.

k:(p,s)∈Wk

s0 = s − (ap su + bp sv + cp , 0)

(8)

as

ρ(s|fp ) = (1 − α) min(kIL (s) − IR (s0 )k1 , τcol )+ α min(|∇x IL (s) − ∇x IR (s0 )|, τgrad ).

(9)

Smoothness term
For the smoothness term, we use a curvature-based, secondorder smoothness regularization term [33] defined as(10)Here, wpq is a contrast-sensitive weight defined as

wpq = e−kIL (p)−IL (q)k1 /γ ,(11)

where γ is a user-defined parameter. The  is a small
constant value that gives a lower bound to the weight ωpq
to increase the robustness to image noise. The functionψ̄pq (fp , fq ) penalizes the discontinuity between fp and fq
in terms of disparity asψ̄pq (fp , fq ) = |dp (fp )−dp (fq )| + |dq (fq )−dq (fp )|,(b) fp = fq

Lemma 1. The term ψpq (fp , fq ) in Eq. (10) satisfies the submodularity of expansion moves in Eq. (2).
Proof. See [33] and also Appendix A.

Here, ∇x I represents the x-component of the gray-value
gradient of image I , and α is a factor that balances theweights of color and gradient terms. The two terms aretruncated by τcol and τgrad to increase the robustness for
occluded regions. We use linear interpolation for IR (s0 ), anda sobel filter kernel of [−0.5 0 0.5] for ∇x . When the data
term is defined on the right image, we swap IL and IR in
Eqs. (7) and (9), and add the disparity value in Eq. (8).

ψpq (fp , fq ) = max (wpq , ) min(ψ̄pq (fp , fq ), τdis ).

(a) fp 6= fq(12)where dp (fq ) = aq pu + bq pv + cq . The first term in Eq. (12)measures the difference between fp and fq by their disparity
values at p, and the second term is defined similarly at q .
We visualize ψ̄pq (fp , fq ) as red arrows in Fig. 2 (a). Theψ̄pq (fp , fq ) is truncated at τdis to allow sharp jumps indisparity at depth edges.

3.3

Local expansion moves

In this section, we describe the fundamental idea of our
method, local expansion moves, as the main contributionof this paper. We first briefly review the original expansion
move algorithm [8], and then describe how we extend it for
efficiently optimizing continuous MRFs.
The expansion move algorithm [8] is a discrete optimization method for pairwise MRFs of Eq. (5), which iteratively
solves a sequence of the following binary labeling problems

f (t+1) = argmin E(f 0 | fp0 ∈ {fp(t) , α})
f0

(13)

for all possible candidate labels ∀ α ∈ L. Here, the binary
variable fp0 for each pixel p is assigned either its current(t)label fp or a candidate label α. If all the pairwise terms in
E(f ) meet the condition of Eq. (2), then the binary energies
E(f 0 ) in Eq. (13) are submodular and this minimization canthus be exactly solved via GC [8] (subproblem optimal).
Here, it is guaranteed that the energy does not increase:E(f (t+1) ) ≤ E(f (t) ). However, the label space L in our
setting is a three dimensional continuous space (a, b, c);therefore, such an exhaustive approach cannot be employed.
Our local expansion moves extend traditional expansion
moves by two ways; localization and spatial propagation. By
localization, we use different candidate labels α dependingon the locations of pixels p in Eq. (13), rather than using thesame α label for all pixels. This is reasonable because the
distributions of disparities should be different from location
to location, therefore the selection of candidate labels αshould be accordingly different. By spatial propagation, we
incorporate label propagation similar to the PatchMatch

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2017

5

Fig. 3. Illustration of the proposed local expansion moves. The local expansion moves consist of many small α-expansions (or local α-expansions),
which are defined using grid structures such shown in the left figure. These local α-expansions are defined at each grid-cell and applied for 3 × 3neighborhood cells (or expansion regions). In the middle part, we illustrate how each of local α-expansions works. (1) The candidate label α (i.e.,α = (a, b, c) representing a disparity plane d = au + bv + c) is produced by randomly choosing and perturbing one of the currently assigned labels
in its center cell. (2) The current labels in the expansion region are updated by α in an energy minimization manner using GC. Consequently, acurrent label in the center cell can be propagated for its surrounding cells. In the right part, local α-expansions are visualized as small patches on
stacked layers with three different sizes of grid structures. As shown here, using local α-expansions we can localize the scopes of label searchingby their locations. Each layer represents a group of mutually-disjoint local α-expansions, which are performed individually in a parallel manner.

inference [1], [2], [5] into GC optimization, and propagate
currently assigned labels to nearby pixels via GC. The
assumption behind this propagation is that, if a good label
is assigned to a pixel, this label is likely a good estimate
for nearby pixels as well. The localization and spatial propagation together make it possible for us to use a powerfulrandomized search scheme, where we no longer need to
produce initial solution proposals as usually done in thefusion based approach [25]. Below we provide the detailed
descriptions of our algorithm.
Local α-expansions for spatial propagation
We first define a grid structure that divides the image
domain Ω into grid regions Cij ⊂ Ω, which are indexedby 2D integer coordinates (i, j) ∈ Z2 . We refer to eachof these grid regions as a cell. We in this paper simply
assume regular square cells, but this could be extended touse superpixels [18], [29], [40]. At a high level, the size ofcells balances between the level of localization and the rangeof spatial propagation. Smaller sizes of cells can achievefiner localization but result in shorter ranges of spatial
propagation. Later in Sec 3.4 we introduce different sizesof multiple grid structures for balancing these two factors
well, but for now let us focus on using one grid structure.
Given a grid structure, we define a local α-expansion
at each cell (i, j), which we specifically denote as αij expansion. We further define two types of regions for each
αij -expansion: its center region Cij and expansion region




[
Rij = Cij ∪
Cmn ,
(14)


(m,n)∈N (i,j)

i.e., 3 × 3 cells consisting of the center region Cij and its
eight neighbor cells.
In the middle part of Fig. 3, we focus on an expansion
region and illustrate how an αij -expansion works. We first
randomly select a pixel r from the center region Cij , andtake its currently assigned label as (a, b, c) = fr . We thenmake a candidate label αij by perturbing this current label
as αij = (a, b, c)+∆. Finally, we update the current labels of
pixels p in the expansion region Rij , by choosing either theircurrent labels fp or the candidate label αij . Here, similarly

Algorithm 1: I TERATIVE αij - EXPANSION
input : current f , target cell (i, j), perturbation size |∆0 |output: updated f (only labels fp at p ∈ Rij are updated)
1 repeat expansion proposer (propagation):2
αij ← fr with randomly chosen r ∈ Cij ;3
f ← argmin E(f 0 | fp0 ∈ {fp , αij }, p ∈ Rij ) ;4 until Kprop times;
5 repeat RANSAC proposer (optional):6
αij ← RANSAC(f, Cij );7f ← argmin E(f 0 | fp0 ∈ {fp , αij }, p ∈ Rij ) ;8 until KRANS times;
9 repeat randomization proposer (refinement):10
αij ← fr with randomly chosen r ∈ Cij ;11
αij ← αij + ∆0 ;12
f ← argmin E(f 0 | fp0 ∈ {fp , αij }, p ∈ Rij ) ;13|∆0 | ← |∆0 |/2 ;0
14 until Krand times (or |∆ | is sufficiently small);to Eq. (13), we update the partial labeling by minimizing
E(f 0 ) with binary variables: fp0 ∈ {fp , αij } for p ∈ Rij , and
fp0 = fp for p ∈/ Rij . Consequently, we obtain an improved
solution as its minimizer with a lower or equal energy.
Notice that making the expansion region Rij larger thanthe label selection region Cij is the key idea for achieving spatial propagation. We can see this since, in an αij expansion without perturbation (∆ = 0), a current label fp
in the center region Cij can be propagated for its nearby
pixels in Rij as the candidate label αij .
We use this αij -expansion iteratively as shown in Algorithm 1. Similarly to the PatchMatch algorithm [1], thisiterative algorithm has propagation (lines 1–4) and refinement (lines 9–14) steps. In the propagation step, we apply
αij -expansions without perturbation to propagate labelsfrom Cij for Rij . In the refinement step, we apply αij expansions with an exponentially decreasing perturbationsize to refine the labels. Here, propagation and refinement
can be seen as heuristic proposers of a candidate label
αij that is expectedly a good estimate for the region Cijor even its neighborhood Rij . Our local expansion movemethod is conceptually orthogonal to any such heuristics
for generating αij . As such an example, we in lines 5–8 use
a RANSAC proposer by following conventional proposal

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2017

generation schemes used in fusion-based methods [33], [45].
Here, we generate αij by fitting a plane to the current
disparity map f in the region Cij using LO-RANSAC [9].
We perform this iterative αij -expansion at every cell (i, j).
This local α-expansion method has the following usefulproperties. Piecewise linearity: It helps to obtain smooth
solutions. In each αij -expansion, multiple pixels in the
expansion region Rij are allowed to move-at-once to thesame candidate label αij at one binary-energy minimization,
which contrasts to BP that updates only one pixel at once.
Since a label represents a disparity plane here, our method
helps to obtain piecewise linear disparity maps and thus
piecewise planar surfaces as we have discussed in Sec 3.1.
Cost filtering acceleration: We can accelerate the computation of matching costs φp (fp ) in Eq. (6) by using costvolume filtering techniques. We discuss this more in Sec 3.5.
Optimality and parallelizability: With our energy formulation, it is guaranteed that each binary-energy minimizationin Algorithm 1 can be optimally solved via GC. In addition,
we can efficiently perform many αij -expansions in a parallelmanner. We discuss these matters in the following sections.
Mutually-disjoint local α-expansions
While the previous section shows how each local αexpansion behaves, here we discuss the scheduling of local
α-expansions. We need a proper scheduling, because local
α-expansions cannot be simultaneously performed due tothe overlapping expansion regions.
To efficiently perform local α-expansions, we divide
them into groups such that the local α-expansions in eachgroup are mutually disjoint. Specifically, we assign each αij expansion a group index k given byk = 4 (j mod 4) + (i mod 4),(15)

and perform the iterative αij -expansions in one group andanother. As illustrated in Fig. 4, this grouping rule picks αij expansions at every four vertical and horizontal cells into
the same group, and it amounts to 16 groups of mutuallydisjoint local α-expansions. We visualize a group of disjoint
local α-expansions as orange regions in Fig. 5, and also asa single layer of stacks in the right part of Fig. 3 with threedifferent grid structures.
Notice that in each group, we leave gaps of one cellwidth between neighboring local α-expansions. These gaps
are to guarantee submodularity and independence for local
α-expansions. By the submodularity, we can show that our
local α-expansions always produce submodular energies
and can thus be optimally solved via GC. Because of thissubmodularity, we can use a standard GC algorithm [7]instead of an expensive QPBO-GC algorithm [22], which
is usually required in the fusion based approach. By the
independence, we can show that the local α-expansions inthe same group do not interfere with each other. Hence, we
can perform them simultaneously in a parallel manner. The
parallelization of GC algorithms are of interest in computer
vision [30], [37]. Our scheme is simple and can directly use
existing GC implementations. The proof of this submodularity and independence is presented in the next section.

6

Algorithm 2: O PTIMIZATION PROCEDURE
Define three levels of grid structures: H = {h1 , h2 , h3 };Initialize the current solution f randomly;
3 Initialize the perturbation size |∆|;
4 repeat5
foreach grid structure h ∈ H do
6
foreach disjoint group k = 0, 1, ..., 15 do
7
foreach cell (i, j) in the group k do [in parallel]8
Do iterative αij expansion (f, (i, j), |∆|);
9
end
10
end
11
end
12
|∆| ← |∆|/2;
13 until convergence;
14 Do post processing;
1

2

Submodularity and independence
To formally address the submodularity and independenceof local α-expansions, we discuss it using the form of
fusion moves [25]. Let us assume a current solution f anda group of mutually-disjoint αij -expansions to be applied.
Simultaneously applying these αij -expansions is equivalent
to the following fusion-energy minimization:f ∗ = argmin E(f 0 | fp0 ∈ {fp , gp }),(16)

where the proposal solution g is set to gp = αij if p belongs
to any of expansion regions Rij in this group; otherwise,p is in gaps so we assign φp (gp ) an infinite unary cost for
forcing fp0 = fp . This proposal disparity map g is visualized
as a single layer of stacks in the right part of Fig. 3. We prove
the following lemmas:
Lemma 2. Submodularity: the binary energies in Eq. (16) aresubmodular, i.e., all the pairwise interactions in Eq. (16) meet thefollowing submodularity of fusion moves [22], [25]:ψpq (gp , gq ) + ψpq (fp , fq ) ≤ ψpq (fp , gq ) + ψpq (gp , fq ). (17)Proof. Omitted.3
Lemma 3. Independence: the assignments to fp0 and fq0 do not
influence each other, if p and q are in different expansion regions.
Proof. The fp0 and fq0 have interactions if and only ifthere exists a chain of valid pairwise interactions C ={ψ(fs0 0 , fs0 1 ), · · · , ψ(fs0 n−1 , fs0 n )} connecting p = s0 and q =sn . But there is no such chain because C inevitably contains
constant or unary terms at a gap (ψcd and ψac in Fig. 5).
3.4

Optimization procedure

Using the local expansion moves shown in the previous
section, we present the overall optimization procedure inthis section and summarize it in Algorithm 2.
This algorithm begins with defining grid structures.
Here, we use thee different sizes of grid structures for betterbalancing localization and spatial propagation.
3. We outline proof. For (p, q) inside an expansion region such as(a, b) in Fig. 5, Eq. (17) is relaxed to Eq. (2). For the other cases, e.g., (a, c)or (c, d), pairwise terms ψac (fa0 , fc0 ) and ψcd (fc0 , fd0 ) become unary or
constant terms. Therefore, when implementing an αij -expansion using
min-cut in Rij , we create nodes for ∀ p ∈ Rij and add ψac (fa0 , fc ) as
unary potentials of nodes a at the inner edge of Rij . See also [8] for theconversion of ψab (fa0 , fb0 ) into edge capacities under expansion moves.

01

Cell index �
2 3 4 5

6

7

0

0

1

2

3

0

1

2

3

0

1

4

5

6

7

4

5

6

7

1

2

8

9

10 11

8

9

10 11

2

3 12 13 14 15 12 13 14 15

0

Cell index �Cell index �IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2017

3

4

0

1

2

3

0

1

2

3

4

5

4

5

6

7

4

5

6

7

6

8

9

10 11

8

9

10 11

6

7 12 13 14 15 12 13 14 15

71

Cell index �
2 3 4 5

6

7

7

8

Matching window ��
Expansion region ����Center region ���

� �
� �

Filtering region ���

�

5

�
�

8

Fig. 4. Group index for αij -expansions.
We perform αij -expansions in the samegroup in parallel.

Fig. 6. Filtering region Mij . The margin width r corresponds with the radius of the matching window Wp .

Fig. 5. Expansion regions of mutually disjoint αij -expansions (group index k = 0).
We leave white gaps between neighbors.

At line 2 of Algorithm 2, the solution f is randomly
initialized. To evenly sample the allowed solution space, wetake the initialization strategy described in [5]. Specifically,for each fp = (ap , bp , cp ) we select a random disparity z0 in
the allowed disparity range [0, dispmax]. Then, a random
unit vector n = (nx , ny , nz ) and z0 are converted to theplane representation by ap = −nx /nz , bp = −ny /nz , and
cp = −(nx pu + ny pv + nz z0 )/nz .
In the main loop through lines 4–13, we select one gridlevel h from the pre-defined grid structures (line 5), andapply the iterative αij expansions of Algorithm 1 for eachcell of the selected structure h (lines 6–10). As discussed inSec. 3.3, we perform the iterative αij expansions by dividing
into disjoint groups defined by Eq. (15). Because of thisgrouping, the αij expansions in the loop at lines 7–9 are
mutually independent and can be performed in parallel.
The perturbation at line 9 of Algorithm 1 is implemented
as described in [5]. Namely, each candidate label αij =(a, b, c) is converted to the form of a disparity d and normal
vector n. We then add a random disparity ∆0d ∈ [−rd , rd ]and a random vector ∆0n of size k∆0n k2 = rn to them,
respectively, to obtain d0 and n0 . Finally, d0 and n0 /|n0 | areconverted to the plane representation αij ← (a0 , b0 , c0 ) toobtain a perturbed candidate label. The values rd and rn
define an allowed change of disparity planes. We initialize
them by setting rd ← dispmax/2 and rn ← 1 at line 3 ofAlgorithm 2, and update them by rd ← rd /2 and rn ← rn /2
at line 12 of Algorithm 2 and line 13 of Algorithm 1.
Finally, we perform post-processing using left-right consistency check and weighted median filtering as describedin [5] for further improving the results. This step is widely
employed in recent methods [3], [5], [14], [32].
Note that there are mainly two differences between thisalgorithm and our previous version [39]. For one thing, wehave removed a per-pixel label refinement step of [39]. This
step is required for updating a special label space structure
named locally shared labels (LSL) used in [39], but can beremoved in our new algorithm by using local α-expansions
instead. For the other, the previous algorithm proceeds
rather in a batch cycle; i.e., it produces all local α-expansions
with all grid structures at once and then applies them tothe current solution f in one iteration. This leads to slower
convergence than our new algorithm that produces each
local α-expansion always from the latest current solution f .

3.5

Fast implementation

Our method has two major computation parts: the calculations of matching costs φp (fp ) of Eq. (6), and application of
GC in Algorithm 1. In Sec. 3.3 and Sec. 3.4, we have shownthat the latter part can be accelerated by performing disjoint
local α expansions in parallel. In this section, we discuss theformer part.
The calculations of the matching costs φp (fp ) are very
expensive, since they require O(|W |) of computation foreach term, where |W | is the size of the matching window. In
our previous algorithm [39], we accelerate this part by using
GPU. The use of GPU is reasonable for our method because
φp (fp ) of all pixels can be individually computed in parallel
during the inference, which contrasts to PMBP [3] that can
only sequentially process each pixel. Still, the computation
complexity is O(|W |) and it becomes very inefficient if weonly use CPU [39].
In the following part, we show that we can approximately achieve O(1) of complexity for computing eachφp (fp ). The key observation here is that, we only needto compute this matching cost φp (fp ) during the αij expansions in Algorithm 1, where φp (fp ) is computed as
φp (αij ) for all p ∈ Rij . With this consistent-label property, we can effectively incorporate the fast subregion costfiltering technique used in [32], by which φp (αij ) for allp ∈ Rij are efficiently computed at once.
To more specifically discuss it, we first separate thecomputation of φp (fp ) into two steps: the calculations of
raw matching costsρfp (s) = ρ(s|fp )(18)

for the support pixels s ∈ Wp , and the aggregation of the
raw matching costs using an edge-aware filter kernel ωpsX
φfp (p) =ωps ρfp (s).
(19)
s∈Wp

We also define a filtering region Mij as the joint matchingwindows in Rij as
[
Mij =
Wp .
(20)
p∈Rij

As shown in Fig. 6, this Mij is typically a square region
margining Rij with r pixels of width around Rij , where r
is the radius of the matching windows.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2017

In the raw matching part of Eq. (18), the calculations of
the raw costs ρfp (s), ∀ s ∈ Wp for all p ∈ Rij generally require O(|W ||Rij |) of total computation. However, with theconsistent-label property (i.e., fp = αij for all p ∈ Rij ), they
can be computed at once in O(|Mij |) by computing ραij (s)for ∀ s ∈ Mij . Here, the computation complexity for each
unary term is O(|Mij |/|Rij |). Therefore, if |Mij | ' |Rij |,we can approximately achieve O(|Mij |/|Rij |) ' O(1) [32].
Similarly, the cost aggregation part of Eq. (19) can bedone in approximately O(1), if we apply a constant-timeedge-aware filtering ωps to ραij (s), ∀ s ∈ Mij [32]. As suchan example we have chosen guided image filtering [13] forEq. (7) but could use more sophisticated filtering [31].
Note that as the consistent-label property is the key tobeing able to use the filtering techniques, this scheme cannot
be used in other PatchMatch based methods [3], [5], [14]except for [32] that uses superpixels as computation units.

48

this calculation by computing each unary term individuallyin parallel on GPU. For both LE-BF and LE-GF, we perform
disjoint local α-expansions in parallel using four CPU cores.
4.1

Evaluation on the Middlebury benchmark V2

We show in Table 1 selected rankings on the Middlebury
stereo benchmark V2 for 0.5-pixel accuracy, comparing withother PatchMatch-based methods [3], [5], [14], [32], [39].
The proposed LE-GF method achieves the current best average rank (3.9) and bad-pixel-rate (5.97%) amongst morethan 150 stereo methods including our previous method(GC+LSL) [39]. Even without post-processing, our LE-GF
method still outperforms the other methods in averagerank, despite that all the methods in Table 1 use thepost-processing. On the other hand, the proposed LE-BF
method achieves comparable accuracy with our previous algorithm [39], since both methods optimize the same energyfunction in very similar ways.

E XPERIMENTS

In the experiments, we first evaluate our method on the
Middlebury benchmark V2 and V3. While the latest V3
benchmark allows us to compare with state-of-the-art methods, we also use the older V2 as there are more related methods [3], [5], [14], [32], [39] registered only in this version.
We further assess the effect of sizes of grid-cells and also
the effectiveness of the proposed acceleration schemes incomparison to our previous method [39]. Our method isfurther compared with the PMBP [3], PMF [32] and Olsson etal. [33] methods that are closely related to our approach. Additional results for these analyses and comparisons on more
image pairs are provided in the supplementary material.
Settings
We use the following settings throughout the experiments.
We use a desktop computer with a Core i7 CPU (3.5 GHz × 4
physical cores) and NVIDIA GeForce GTX Titan Black GPU.
All methods are implemented using C++ and OpenCV.
The parameters of our data term are set as{e, τcol , τgrad , α} = {0.012 , 10, 2, 0.9} following [32], [36] and[5]. The size of matching windows Wp is set to 41 × 41 (i.e.,the size of local regression windows Wk0 in Eq. (7) is set
to 21 × 21), which is the same setting with PMBP [3] and[5], [14]. For the smoothness term, we use {λ, τdis , , γ} ={1, 1, 0.01, 10} and eight neighbors for N .
For our algorithm, we use three grid structures with cell
sizes of 5 × 5, 15 × 15, and 25 × 25 pixels. The iteration
numbers {Kprop , Krand } in Algorithm 1 are set to {1, 7} forthe first grid structure, and {2, 0} (only propagation step)for the other two. Except in Sec. 4.2 and 4.7 we disable
the RANSAC proposer for fair comparisons with related
methods or for avoiding cluttered analysis. When enabled,
we set KRANS = 1 for all the grid structures. We iterate themain loop ten times. We use a GC implementation of [7].
We use two variants of our method. LE-GF uses guided
image filtering [13] for wps as proposed in Sec. 3.2. We only
use a CPU implementation for this method. LE-BF uses
bilateral filtering for ωps defined similarly to Eq. (11) and
λ = 20 so that the energy function corresponds to one weused in [39]. The computation of matching costs is thus asslow as O(|W |). In a GPU implementation, we accelerate

4.2

Evaluation on the Middlebury benchmark V3

As major differences from the older version, the latest Middlebury benchmark introduces more challenging difficulties such as high-resolution images, different exposure andlighting between an image pair, imperfect rectification, etc.
Since our model, especially our data term borrowed from[5], does not consider such new difficulties, we here slightly
modify our model for adapting it to the latest benchmark.
To this end, we incorporate the state-of-the-art CNN-basedmatching cost function by Zbontar and LeCun [49], byfollowing manners of current top methods [10], [20], [28],[29] on the latest benchmark. Here, we only replace our
pixelwise raw matching cost function ρ(·) of Eq. (9) with
their matching cost function as follows.

ρ0 (s|fp ) = min(CCNN (s, s0 ), τCNN ).

(21)

0

The function CCNN (s, s ) computes a matching cost between
left and right image patches of 11×11-size centered at s and
s0 , respectively, using a neural network (called MC-CNNacrt in [49]).4 In our method we truncate the matching costvalues at τCNN = 0.5 and aggregate them for the support
window pixels s ∈ Wp using slanted patch matching of
Eq. (6). In addition to the change in Eq. (9), the widths
of square cells of three grid structures are changed proportionally to 1%, 3%, and 9% of the image width. Since
only one submission is allowed for the test data by the
new benchmark regulation, we here only use our LE-GF
method with λ = 0.5 and KRANS = 1. We keep all the otherparameters (i.e., {e, Wp , τdis , , γ, Kprop , Krand }) as default.
Table 2 shows selected rankings for the bad 2.0 metric(i.e., percentage of bad pixels by the error threshold of
2.0 pixels at full resolution) for non-occluded regions as
the default metric. Our method outperforms all existing 63
methods on the table not only for the default metric but
for all combinations of {bad 0.5, 1.0, 2.0, 4.0} × {nonocc, all}
4. Using the authors’ code and pre-trained model, we pre-computethe matching costs CCNN at all allowed integer disparities using frontoparallel patches. During the inference, we linearly interpolate thefunction CCNN (s, s0 ) at non-integer coordinates s0 . MC-CNN learns
a patch similarity metric invariantly to the left-right patch distortion.
Thus, the bias due to using fronto-parallel patches in CCNN is small.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2017

9

TABLE 1
Middlebury benchmark V2 for 0.5-pixel accuracy. Our method using the guided image filtering (GF) achieves the current best average rank 3.9.
Error percents for all pixels, non-occluded pixels (nonocc), and pixels around depth discontinuities (dis) are shown. (Snapshot on April 22, 2015)

TABLE 2
Middlebury benchmark V3 for the bad 2.0 nonocc metric. Our method achieves the current best average error rate among 64 existing algorithms.
We here list current top methods [10], [20], [28], [29], [49] that use the matching costs of MC-CNN-acrt [49] like ours. (Snapshot on July 4, 2017)

except for bad 4.0 - all. In Table 3 we analyze effects ofour post-processing and RANSAC proposer using training
data. Again, our method is ranked first even without postprocessing for the default and all other bad nonocc metrics.
Also, the RANSAC proposer reduces errors by one point(see Sec. 4.7 for more analysis). The avr and stdev show that
our inference is stable by different random initializations.
4.3

Effect of grid-cell sizes

To observe the effect of grid-cell sizes, we use the three sizesof grid-cells, 5×5 pixels (denoted as “S”mall), 15×15 pixels(denoted as “M”edium), 25×25 pixels (denoted as “L”arge),
in different combinations and assess the performance using
the following five different settings; (S, S, S): the small-size
cells for all grid-structures; (M, M, M): the medium-size cells
for all grid-structures; (L, L, L): the large-size cells for all
grid-structures; (S, M, M): the small and medium-size cells
for the first and the other two grid-structures, respectively;(S, M, L): the small, medium, and large-size cells for the first,
second, and third grid-structures, respectively (the default
setting for our method described above). Here, the iteration
numbers for the first to third grid-structures are kept as
default. We use the LE-GF method so as to access the effecton cost filtering acceleration as well. We use λ = 0.5 but
keep the other parameters as default. Using these settings,
we observe the performance variations by estimating the
disparities of only the left image of the Reindeer dataset
without post-processing.
The plots in Fig. 7 (a) show the transitions of the energy function values over iterations. The plots in Fig. 7 (b)show the temporal transitions of error rates with subpixel
accuracy. As shown, the joint use of multiple cell sizes
improves the performance in both energy reduction andaccuracy. Although (S, M, M) and (S, M, L) show almost
the same energy transitions, the proposed combination (S,M, L) shows faster and better convergence in accuracy.
Figure 8 compares the results of the five settings witherror maps for (S, M, M) and (S, M, L). The use of larger gridcells helps to obtain smoother disparities, and it is especiallyeffective for occluded regions.
Comparing the running times in Fig. 7 (b), the use of thesmall-size cells is inefficient due to the increased overheadin cost filtering, whereas the use of the medium and largesize cells achieves almost the same efficiency.
4.4 Efficiency evaluation in comparison to our previous algorithm [39]We evaluate the effectiveness of the three acceleration techniques: 1) parallelization of disjoint α-expansions, and 2)
acceleration of unary cost computation by GPU and 3) bycost filtering. For this, we compare following six variants ofour method: LE-BF and LE-GF using one or four CPU cores(denoted as CPUx1 and CPUx4), and LE-BF using GPU
and one or four CPU cores (denoted as GPU+CPUx1 andGPU+CPUx4). Additionally, we compare with our previous
algorithm of [39] with CPU and GPU implementations(denoted as LSL with CPUx1 and GPU+CPUx4). We use
the Rocks1 dataset by estimating the disparities of only theleft image without post-processing. Figures 9 (a)–(c) show

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2017
0.57

ML,bL,bLL

0.54
0.53

avr ± stdev

X

X

ref. 3DMST [28]

X

Ours

nonocc
6.52
6.65
7.72
6.63
±0.12
7.08

all
12.1
13.6
14.6
12.3
±0.2
12.9

23

MS,bM,bML

MS, S,bSL

MS,bM,bLL

0.52
MM, M,bML

0.51
0.5

ML, L,bLL

0.48(S, S,,S)

21(L, L,,L)20(M, M,,M)(S, M,,M)(S, M,,L)18

MS, M,bML
MS,bM,bLL
1

22

19

0.49

0

RP
X
X

24

MM,bM,bML

0.55

0.47

PP
X

25

MS,bS,bSL

0.56

Error,Rate,[%]RelativebEnergybFunctionbValueTABLE 3
Effect of our post-processing (PP) and RANSAC
proposer (RP). We show weighted average bad
2.0 scores of our LE-GF method with 3DMST [28]that ranks second following ours. We use 15
training image pairs from Middlebury V3. We also
show average (avr) and standard deviation(stdev) by ten sets of different random
initializations. See also Sec. 4.7 for more analysis
on the RANSAC proposer.

10

17

2

3

4

5

6

7

8

9

10

0

50

NumberbofbIterations

100

150

200

Running,Time,[seconds](a) Energy function values w.r.t. iterations(b) Temporal error-rate transitions

Fig. 7. Effect of grid-cell sizes. We use LE-GF with different combinations of grid structures. The
S, M, and L denote small, medium, and large grid-cells, respectively. The joint use of different
sizes of grid-cells improves the performance. See also Fig. 8 for visual comparison.

Reindeer(S, S, S)(S, M, M)(S, M, L)

Ground truth(M, M, M)(L, L, L)

Error map of (S, M, M)

Error map of (S, M, L)

Visibility map

Fig. 8. Visual effect of grid-cell sizes. The use of larger grid-cells leads to smoother solutions and effective for occluded regions. The proposed
combination (S, M, L) well balances localization and spatial propagation and performs best. These are all raw results without post-processing.
0+68

LSLsPCPUx1U
Puntils6316ssecU

0.90
0.85
0.80

0.75
0.70

LE3BFsPCPUx1U
Puntils1894ssecU

0.65
0.60

28

LSLsCGPUxCPUx4P

24
0+66
0+65

LESGF CCPUx1P
CCPUx4P

0+64

LESBF CCPUx1P
CCPUx4P

200

400

600

800

RunningsTimes[seconds](a) Time-energy transitions (full)

1000

22
LE-BF
FCPUx4)20

LE-BFPFGPU+CPUx1)
FGPU+CPUx4)

14

LE-GF
FCPUx4)P FCPUx1)

12
0

100

200

300

400

RunningsTimes[seconds](b) Time-energy transitions (zoomed)

LSLPFCPUx1)
LSLPFGPU+CPUx4)
LE-BFPFCPUx1)
LE-BFPFCPUx4)
LE-GFPFCPUx1)
LE-GFPFCPUx4)
LE-BLPFGPU+CPUx1)
LE-BLPFGPU+CPUx4)18
16

LESGF CGPUxCPUx1P
CGPUxCPUx4P

0+63
0+62

0

LSL
FGPU+CPUx4)

26

0+67

ErrorPRateP[+]

0.95

RelativesEnergysFunctionsValue

RelativesEnergysFunctionsValue

1.00

0

100

200

300

400

RunningPTimeP[seconds](c) Error rate transitions

Fig. 9. Efficiency evaluation in comparison to our previous algorithm (LSL) [39]. Accuracies are evaluated for all-regions at each iteration.

the temporal transitions of the energy values in the full andzoomed scales, and the subpixel error rates, respectively.
Note that energy values in Fig. 9 are evaluated by the energyfunction of LE-BF and LSL for better visualization.
Parallelization of local α-expansions on CPU: Comparing CPUx1 and CPUx4, we observe about 3.5x of speed-up
for both LE-BF and LE-GF. Note that the converged energy
values of LE-GF are relatively higher than the others because
it optimizes a different energy function. However, if we see
Fig. 9 (c), it actually finds better solutions in this example.

On the other hand, speed-up for LE-BF from GPU+CPUx1
to GPU+CPUx4 is limited to about 1.7x. This is because theunary cost computation is already parallelized by GPU andthis speed-up is accounted for by the parallelization of otherparts, e.g., the computation of pairwise terms and min-cut.
Parallelization of unary cost computation on GPU: By
comparing GPU+CPUx1 and CPUx1 of LE-BF, we obtainabout 19x of speed-up. This is relatively smaller than 32x ofspeed-up in our previous LSL method, mainly due to thelarger overhead of GPU data transferring.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2017

Fast cost filtering: By comparing LE-GF and LE-BF
methods, we observe about 5.3x speed-up for both CPUx1
and CPUx4 versions. Notably, even a CPU implementationof LE-GF (CPUx4) achieves almost comparable efficiency
with a GPU implementation of LE-BF (GPU+CPUx1).
Comparison to our previous method [39]: In comparison to our previous LSL method [39], our LE-BF showsmuch faster convergence than LSL in the same singlecore results (CPUx1). We consider there are mainly threefactors contributing to this speed-up; First, the batch-cycle
algorithm adopted in [39] makes its convergence slower;
Second, we use cost filtering and compute its first stepof Eq. (18) in O(1), while [39] computes each unary term
individually in O(|W |); Finally, we have removed a perpixel label refinement step of [39] which requires additional
unary-cost computations. Similar speed-up can be observed
from LSL (GPU+CPUx4) to LE-BF (GPU+CPUx4).

pixels K is set to 300, 500, and 700 as used in [32] (K = 500
is the default in [32]), and we sufficiently iterate 30 times.
Both LE-GF and PMF are run using a single CPU core.
Figures 12 (a)–(c) show the temporal transitions of theenergy values in the full and zoomed scales, and the subpixel error rates, respectively. As shown in Figs. 12 (a)and (b), PMF converges at higher energies than ours, since
it cannot explicitly optimize pairwise smoothness terms
because it is a local method. Furthermore, although energies
are reduced almost monotonically in PMF, the transitions ofaccuracies are not stable and even degrade after many iterations. This is also because of the lack of explicit smoothness
regularizer, and PMF converges at a bad local minimum.
Figure 13 compares the results of our method and PMF with
K = 500. Again, our methods find smoother and better
disparities around edges and at occluded regions.
4.7

4.5

Comparison with PMBP [3]We compare our method with PMBP [3] that is the closest
method to ours. For a fair comparison, we use four neighbors for N in Eq. (5) as the same setting with PMBP. For
a comparable smoothness weight with the default setting(eight-neighbor N ), we use λ = 40 for LE-BF and λ = 2
for LE-GF, and keep the other parameters as default. For
PMBP, we use the same model as ours; the only difference
from the original PMBP is the smoothness term, which
does not satisfy the submodularity of Eq. (2). In PMBP, it
defines K candidate labels for each pixel, for which we set
K = 1 and K = 5 (original paper uses K = 5). We show
the comparison using the Cones dataset by estimating the
disparities of only the left image without post-processing.
Figures 10 (a)–(c) show the temporal transitions of theenergy values in the full and zoomed scales, and the subpixel error rates, respectively. We show the performanceof our method using its GPU and CPU (one or four CPU
cores) implementations. For PMBP, we also implemented
the unary cost computation on GPU, but it became rather
slow, due to the overhead of data transfer. Efficient GPU
implementations for PMBP are not available in literature.5
Therefore, the plots show PMBP results that use a single
CPU core. Figures 10 (a) and (b) show that, even with a single CPU-core implementation, our LE-BF and LE-GF show
comparable or even faster convergence than PMBP. WithCPU and GPU parallelization, our methods achieve much
faster convergence than PMBP. Furthermore, our methods
reach lower energies with greater accuracies than PMBP at
the convergence. In Fig. 11 we compare the results of our
LE-GF method and PMBP with K = 5. While PMBP yields
noisy disparities, our method finds smoother and better
disparities around edges and at occluded regions.
4.6

Comparison with PMF [32]We also compare our LE-GF method with PMF [32] using
the same data term as ours. For PMF, the number of super5. GPU-parallelization schemes of BP are not directly applicable due
to PMBP’s unique settings. The “jump flooding” used in the original
PatchMatch [1] reports 7x speed-ups by GPU. However, because itpropagates candidate labels to distant pixels, it is not applicable to
PMBP that must propagate messages to neighbors, and is not as efficientas our 19x, either.

11

Comparison with Olsson et al. [33]

We compare with a stereo method by Olsson et al. [33], fromwhich we borrow our smoothness term. As differences fromour method, their data term does not use accurate slanted
patch matching and they use conventional fusion-based optimization using planar proposals generated by RANSACbased local plane fitting (see [33] for details). To compare
optimization approaches of ours and [33], we replace thedata term of [33] with ours so that the two methods use
the same energy function as ours. We here use two variantsof our LE-GF method that remove the RANSAC proposer
and the regularization term, which are components closely
related to [33]. We disable post-processing in all methods.
Table 4 shows error rates for Teddy from Middlebury
V2 and three more datasets from V3, where the proposed
method outperforms [33] and others. As shown in Figs. 14
and 15, our method has faster convergence and produces
qualitatively better smooth surfaces than [33]. Spatial label
propagation not only accelerates inference but also introduces an implicit smoothness effect [5], which helps our
method to find smoother surfaces. Our method mostly
performs comparably well regardless of the RANSAC proposer, outperforming [33]. On the other hand, the RANSAC
proposer promotes convergence (see Fig. 14) and also effectively recovers large texture-less planes in Vintage (Fig. 16).
Using a more advanced proposer [34] may lead to further
improvements. Also note that, except for Vintage in Table 4,
our method performs well for non-occluded regions even
without regularization (as comparable to [33]). The slanted
patch matching term alone has an implicit second-ordersmoothness effect [5], and the regularizer of [33] further
enforces it especially for occluded or texture-less regions.

5

C ONCLUSIONSIn this paper, we have presented an accurate and efficientstereo matching method for continuous 3D plane label
estimation. Unlike previous approaches that use fusion
moves [25], [33], [45], our method is subproblem optimaland only requires a randomized initial solution. By comparing with a recent continuous MRF stereo method, PMBP [3],our method has shown an advantage in efficiency andcomparable or greater accuracy. The use of a GC-based
optimizer makes our method advantageous.

1x00

GK8+

0x95

GK79

0x90
0x85

PMBP =K=5)0x80

LEUGF =CPUx1)

0x75

LEUBF =CPUx1)
PMBP =K=1)

0x70
0x65

12
20

PMBP[CK=+P

18

GK77
GK75

LEMGF
CCPUx+P
CCPUx4P

GK73
GK7+

LEMBF
CGPUxCPUx4P

GK69

LEMBF CCPUx+P

300

600

900

1200

1500

17
16
15

PMBPMFK=5C
F12.4)MafterM2MiterC

14

LE+BF
FGPUxCPUx4C

13
12
10

G

RunningsTimes[seconds]

+GG

2GG

3GG

4GG

5GG

0

Running[Time[[seconds](a) Time-energy transitions (full)

PMBPMFK=1C

LE+GF
FCPUx1CFCPUx4C

11

GK67
0

PMBPMFK=1C
PMBPMFK=5C
LE+BFMFCPUx1C
LE+GFMFCPUx1C
LE+GFMFCPUx4C
LE+BLMFGPUxCPUx4C

19

ErrorMRateM[)]

Relative[Energy[Function[Value

RelativesEnergysFunctionsValue

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2017300

LE+BF FCPUx1C
600

900

1200

1500

RunningMTimeM[seconds](b) Time-energy transitions (zoomed)(c) Error rate transitions

Fig. 10. Efficiency and accuracy comparison with PMBP [3]. Our methods achieve much faster convergence, reaching lower energies and better
accuracies at the convergence. Accuracies are evaluated for all-regions at each iteration. See also Fig. 11 for visual comparison.

(a) Our LE-BF(b) PMBP [3](c) Ground truth

Fig. 11. Visual comparison with PMBP [3]. Our method finds smother and better disparities around edges and at occluded regions.
0.57

0.90
0.80

0.70
0.60

PMF-GF

LE-GF

0.50

15

PMF-GF

0.55

0.54
0.53

100

200

300

Running[Time[[seconds](a) Time-energy transitions (full)

400

13
12

PMF-GFGUK=700x
PMF-GFGUK=500x
PMF-GFGUK=300x
LE-GFGUCPUx1x

PMF-GF
UK=700x

PMF-GF
UK=300x

11

10

LE-GF

0.52
0

PMF-GF
UK=500x140.56

ErrorGRateG[)]

RelativedEnergydFunctiondValue

Relative[Energy[Function[Value

1.00

LE-GF

9
0

100

200

300

400

RunningdTimed[seconds](b) Time-energy transitions (zoomed)

0

100

200

300

400

RunningGTimeG[seconds](c) Error rate transitions

Fig. 12. Efficiency and accuracy comparison with PMF [32]. Our method stably improves the solution and reaches a lower energy with greater
accuracy at the convergence. Accuracies are evaluated for all-regions at each iteration. See also Fig. 13 for visual comparison.

(a) Our LE-GF(b) PMF [32](c) Ground truth

Fig. 13. Visual comparison with PMF [32]. With explicit smoothness regularization, our method finds smoother and better disparities around edges
and at occluded regions.

Furthermore, by using the subregion cost-filtering
scheme developed in a local stereo method PMF [32], we
achieve a fast algorithm and greater accuracy than PMF. As
shown in [32], our method will be even faster by using amore light-weight constant-time filtering such as [31].

Follow-up works [18], [40] suggest that our optimization
strategy can be used for more general correspondence field
estimation than stereo. We also believe that occlusion handling schemes using GC [23], [44] can be incorporated into
our framework, which may yield even greater accuracy.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2017

13

TABLE 4
Accuracy comparison with Olsson et al. [33]. We use our energy function using 3D patch matching for our method and [33]. Teddy dataset is from
Middlebury V2 (Sec. 4.1) and the other three from V3 (Sec. 4.2). Our local expansion move method (LE) outperforms fusion-based optimization of[33]. While our method has comparably good accuracies regardless of the RANSAC proposer, it is effective for large texture-less planes in Vintage.

Ours (proposed LE-GF)
Ours w/o RANSAC
Ours w/o regularization
Olsson et al. [33] + 3DPM

Reg.
[33]
[33]
None
[33]

Options
Optimization
LE+RANSACLE
LE+RANSAC
Fusion+RANSAC

Teddy (from V2)
nonocc
all
disc
3.98
9.46
12.8
4.26
9.19
13.6
5.47
13.0
16.0
5.21
9.98
15.8

Adirondack
nonocc
all
1.19
4.36
1.20
4.63
2.56
7.74
2.74
5.42

ArtL
nonocc
all
3.97
15.5
3.98
15.3
4.12
17.6
4.41
15.8

Vintagenonocc
all
5.65
11.7
13.7
18.7
22.8
27.9
5.94
11.8

Error Rate [%]18
16

Ours (proposed)

14

Ours w/o RANSAC

12

Olsson et al.

10

Olsson et al.

8
6
4

Ours
(proposed)

2
0
0

500

1000

1500

VintageOurs (proposed)

Ours w/o RANSAC

Olsson et al. [33]Ours w/o
RANSAC
2000

2500

3000

3500

Running Time [seconds](For Andirondack, by a single CPU core)Fig. 14. Efficiency and accuracy comparison
with Olsson et al. [33]. Our method is fasterthan [33] even without parallelization. TheRANSAC proposer promotes the inference.

Our LE-GF

Fig. 15. 3D visual comparison with Olsson etal. [33] for ArtL dataset. Although we use thesame energy function, the result by [33] isstrongly biased to piecewise planar surfaces.

ACKNOWLEDGMENTS
The authors would like to thank the anonymous reviewers
including those of the conference paper. This work was
supported by JSPS KAKENHI Grant Number 14J09001 andMicrosoft Research Asia Ph.D. Fellowship.

R EFERENCES[1]

Olsson et al. [33]

C. Barnes, E. Shechtman, A. Finkelstein, and D. B. Goldman.
PatchMatch: a randomized correspondence algorithm for structural image editing. Proc. of SIGGRAPH (ACM Trans. on Graph.),28(3):24:1–24:11, 2009.
[2] C. Barnes, E. Shechtman, D. B. Goldman, and A. Finkelstein. The
generalized patchmatch correspondence algorithm. In Proc. ofEuropean Conf. on Computer Vision (ECCV), pages 29–43, 2010.
[3] F. Besse, C. Rother, A. Fitzgibbon, and J. Kautz. PMBP: PatchMatch
Belief Propagation for Correspondence Field Estimation. In Proc.
of British Machine Vision Conf. (BMVC), pages 132.1–132.11, 2012.
[4] S. Birchfield and C. Tomasi. Multiway cut for stereo and motion
with slanted surfaces. In Proc. of Int’l Conf. on Computer Vision(ICCV), volume 1, pages 489–495, 1999.
[5] M. Bleyer, C. Rhemann, and C. Rother. PatchMatch Stereo Stereo Matching with Slanted Support Windows. In Proc. of British
Machine Vision Conf. (BMVC), pages 14.1–14.11, 2011.
[6] M. Bleyer, C. Rother, and P. Kohli. Surface stereo with soft
segmentation. In Proc. of IEEE Conf. on Computer Vision and Pattern
Recognition (CVPR), pages 1570–1577, 2010.
[7] Y. Boykov and V. Kolmogorov. An experimental comparison ofmin-cut/max- flow algorithms for energy minimization in vision.
IEEE Trans. Pattern Anal. Mach. Intell., 26(9):1124–1137, 2004.
[8] Y. Boykov, O. Veksler, and R. Zabih. Fast approximate energy
minimization via graph cuts. IEEE Trans. Pattern Anal. Mach. Intell.,23(11):1222–1239, 2001.
[9] O. Chum, J. Matas, and J. Kittler. Locally Optimized RANSAC. In
DAGM-Symposium, volume 2781, pages 236–243, 2003.
[10] S. Drouyer, S. Beucher, M. Bilodeau, M. Moreaud, and L. Sorbier.
Sparse Stereo Disparity Map Densification using Hierarchical Image Segmentation. In Proc. of Int’l Symp. on Mathematical Morphology, pages 172–184, 2017.

Fig. 16. Effectiveness of the RANSAC-basedplane proposer for large texture-less regions.
We show error maps of our method and Olsson et al. [33] for Vintage dataset.

[11] P. Felzenszwalb and D. Huttenlocher. Efficient belief propagationfor early vision. In Proc. of IEEE Conf. on Computer Vision and
Pattern Recognition (CVPR), volume 1, pages 261–268, 2004.
[12] S. Geman and D. Geman. Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images. IEEE Trans. PatternAnal. Mach. Intell., 6(6):721–741, 1984.
[13] K. He, J. Sun, and X. Tang. Guided image filtering. IEEE Trans.
Pattern Anal. Mach. Intell., 35(6):1397–1409, 2013.
[14] P. Heise, S. Klose, B. Jensen, and A. Knoll. PM-Huber: PatchMatchwith Huber Regularization for Stereo Matching. In Proc. of Int’l
Conf. on Computer Vision (ICCV), pages 2360–2367, 2013.
[15] L. Hong and G. Chen. Segment-based stereo matching usinggraph cuts. In Proc. of IEEE Conf. on Computer Vision and Pattern
Recognition (CVPR), volume 1, pages 74–81, 2004.
[16] A. Hosni, M. Gelautz, and M. Bleyer. Accuracy-Efficiency Evaluation of Adaptive Support Weight Techniques for Local Stereo
Matching. In Proc. DAGM/OAGM Symposium, pages 337–346, 2012.
[17] A. Hosni, C. Rhemann, M. Bleyer, C. Rother, and M. Gelautz. Fast
cost-volume filtering for visual correspondence and beyond. IEEE
Trans. Pattern Anal. Mach. Intell., 35(2):504–511, 2013.
[18] J. Hur and S. Roth. MirrorFlow: Exploiting Symmetries in Joint
Optical Flow and Occlusion Estimation. In Proc. of Int’l Conf. on
Computer Vision (ICCV), 2017. (to appear).
[19] H. Ishikawa. Exact optimization for Markov random fields withconvex priors. IEEE Trans. Pattern Anal. Mach. Intell., 25(10):1333–
1336, 2003.
[20] K. R. Kim and C. S. Kim. Adaptive smoothness constraints for
efficient stereo matching using texture and edge information. In
Proc. of Int’l Conf. on Image Processing (ICIP), pages 3429–3433, 2016.
[21] V. Kolmogorov. Convergent Tree-Reweighted Message Passing
for Energy Minimization. IEEE Trans. Pattern Anal. Mach. Intell.,28(10):1568–1583, 2006.
[22] V. Kolmogorov and C. Rother. Minimizing Nonsubmodular Functions with Graph Cuts – A Review. IEEE Trans. Pattern Anal. Mach.
Intell., 29(7):1274–1279, 2007.
[23] V. Kolmogorov and R. Zabih. Computing visual correspondencewith occlusions using graph cuts. In Proc. of Int’l Conf. on Computer
Vision (ICCV), volume 2, pages 508–515, 2001.
[24] V. Kolmogorov and R. Zabin. What energy functions can be
minimized via graph cuts? IEEE Trans. Pattern Anal. Mach. Intell.,26(2):147–159, 2004.
[25] V. Lempitsky, C. Rother, S. Roth, and A. Blake. Fusion Moves

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2017

for Markov Random Field Optimization. IEEE Trans. Pattern Anal.
Mach. Intell., 32(8):1392–1405, 2010.
[26] A. Li, D. Chen, Y. Liu, and Z. Yuan. Coordinating multiple
disparity proposals for stereo computation. In Proc. of IEEE Conf.
on Computer Vision and Pattern Recognition (CVPR), pages 4022–
4030, 2016.
[27] G. Li and S. W. Zucker. Differential geometric inference in surface
stereo. IEEE Trans. Pattern Anal. Mach. Intell., 32(1):72–86, 2010.
[28] L. Li, X. Yu, S. Zhang, X. Zhao, and L. Zhang. 3d cost aggregation
with multiple minimum spanning trees for stereo matching. Appl.
Opt., 56(12):3411–3420, 2017.
[29] L. Li, S. Zhang, X. Yu, and L. Zhang. PMSC: PatchMatch-Based
Superpixel Cut for Accurate Stereo Matching. IEEE Transactions on
Circuits and Systems for Video Technology, PP(99):1–1, 2017.
[30] J. Liu and J. Sun. Parallel graph-cuts by adaptive bottom-up
merging. In Proc. of IEEE Conf. on Computer Vision and Pattern
Recognition (CVPR), pages 2181–2188, 2010.
[31] J. Lu, K. Shi, D. Min, L. Lin, and M. Do. Cross-based local
multipoint filtering. In Proc. of IEEE Conf. on Computer Vision and
Pattern Recognition (CVPR), pages 430–437, 2012.
[32] J. Lu, H. Yang, D. Min, and M. N. Do. Patch Match Filter: Efficient
Edge-Aware Filtering Meets Randomized Search for Fast Correspondence Field Estimation. In Proc. of IEEE Conf. on Computer
Vision and Pattern Recognition (CVPR), pages 1854–1861, 2013.
[33] C. Olsson, J. Ulen, and Y. Boykov. In Defense of 3D-Label Stereo.
In Proc. of IEEE Conf. on Computer Vision and Pattern Recognition(CVPR), pages 1730–1737, 2013.
[34] C. Olsson, J. Uln, and A. Eriksson. Local Refinement for Stereo
Regularization. In Proc. of Int’l Conf. on Pattern Recognition (ICPR),
pages 4056–4061, 2014.
[35] G. Petschnigg, R. Szeliski, M. Agrawala, M. Cohen, H. Hoppe, andK. Toyama. Digital photography with flash and no-flash image
pairs. ACM Trans. on Graph., 23(3):664–672, 2004.
[36] C. Rhemann, A. Hosni, M. Bleyer, C. Rother, and M. Gelautz.
Fast cost-volume filtering for visual correspondence and beyond.
In Proc. of IEEE Conf. on Computer Vision and Pattern Recognition(CVPR), pages 3017–3024, 2011.
[37] P. Strandmark and F. Kahl. Parallel and distributed graph cuts by
dual decomposition. In Proc. of IEEE Conf. on Computer Vision and
Pattern Recognition (CVPR), pages 2085–2092, 2010.
[38] R. Szeliski, R. Zabih, D. Scharstein, O. Veksler, V. Kolmogorov,
A. Agarwala, M. Tappen, and C. Rother. A Comparative Studyof Energy Minimization Methods for Markov Random Fields with
Smoothness-Based Priors. IEEE Trans. Pattern Anal. Mach. Intell.,30(6):1068–1080, 2008.
[39] T. Taniai, Y. Matsushita, and T. Naemura. Graph Cut based
Continuous Stereo Matching using Locally Shared Labels. In Proc.
of IEEE Conf. on Computer Vision and Pattern Recognition (CVPR),
pages 1613–1620, 2014.
[40] T. Taniai, S. N. Sinha, and Y. Sato. Joint Recovery of Dense
Correspondence and Cosegmentation in Two Images. In Proc.
of IEEE Conf. on Computer Vision and Pattern Recognition (CVPR),
pages 4246–4255, 2016.
[41] H. Tao, H. Sawhney, and R. Kumar. A global matching frameworkfor stereo computation. In Proc. of Int’l Conf. on Computer Vision(ICCV), volume 1, pages 532–539, 2001.
[42] L. Wang and R. Yang. Global stereo matching leveraged by sparse
ground control points. In Proc. of IEEE Conf. on Computer Vision
and Pattern Recognition (CVPR), pages 3033–3040, 2011.
[43] Z.-F. Wang and Z.-G. Zheng. A region based stereo matchingalgorithm using cooperative optimization. In Proc. of IEEE Conf.
on Computer Vision and Pattern Recognition (CVPR), 2008.
[44] Y. Wei and L. Quan. Asymmetrical Occlusion Handling Using
Graph Cut for Multi-View Stereo. In Proc. of IEEE Conf. on
Computer Vision and Pattern Recognition (CVPR), pages 902–909,
2005.
[45] O. Woodford, P. Torr, I. Reid, and A. Fitzgibbon. Global Stereo
Reconstruction under Second-Order Smoothness Priors. IEEE
Trans. Pattern Anal. Mach. Intell., 31(12):2115–2128, 2009.
[46] K. Yamaguchi, T. Hazan, D. McAllester, and R. Urtasun. Continuous markov random fields for robust stereo estimation. In Proc. ofEuropean Conf. on Computer Vision (ECCV), pages 45–58, 2012.
[47] J. S. Yedidia, W. T. Freeman, and Y. Weiss. Generalized Belief
Propagation. In Advances in Neural Information Processing Systems(NIPS), volume 13, pages 689–695, 2000.
[48] K.-J. Yoon and I.-S. Kweon. Locally adaptive support-weightapproach for visual correspondence search. In Proc. of IEEE Conf.
on Computer Vision and Pattern Recognition (CVPR), volume 2, pages
924–931, 2005.
[49] J. Zbontar and Y. LeCun. Stereo matching by training a con-14

volutional neural network to compare image patches. Journal of
Machine Learning Research, 17:1–32, 2016.

Tatsunori Taniai received his B.E., M.S., andPh.D. degrees in EECS from the University ofTokyo in 2012, 2014, 2017, respectively. Since
2017 he is with RIKEN Center for AdvancedIntelligence Project (AIP) in Japan as a postdoctoral researcher. His research interests lie in low
and mid-level computer vision including image
based 3D scene reconstruction, Markov random
field optimization, and image segmentation.

Yasuyuki Matsushita received his B.S., M.S.
and Ph.D. degrees in EECS from the Universityof Tokyo in 1998, 2000, and 2003, respectively.
From April 2003 to March 2015, he was with
Visual Computing group at Microsoft Research
Asia. In April 2015, he joined Osaka University
as a professor. His research area includes computer vision, machine learning and optimization.
He is on the editorial board of IEEE Transactions
on Pattern Analysis and Machine Intelligence(TPAMI), International Journal of Computer Vision (IJCV), IPSJ Journal of Computer Vision and Applications (CVA),The Visual Computer Journal, and Encyclopedia of Computer Vision.
He served/is serving as a Program Co-Chair of PSIVT 2010, 3DIMPVT
2011, ACCV 2012, ICCV 2017, and a General Co-Chair for ACCV 2014.
He is a senior member of IEEE.

Yoichi Sato received the BS degree from the
University of Tokyo in 1990, and the MS and PhD
degrees in robotics from the School of Computer
Science, Carnegie Mellon University, in 1993
and 1997, respectively. He is a professor at theInstitute of Industrial Science, the University ofTokyo. His research interests include physicsbased vision, reflectance analysis, and gaze andgesture analysis. He served/is serving to conference organization and journal editorial roles
including IEEE Transactions on Pattern Analysis
and Machine Intelligence, International Journal of Computer Vision,
IPSJ Journal of Computer Vision and Applications, ECCV 2012 program
co-chair, MVA 2013 general chair, ACCV 2018 general co-chair, and3DV 2018 program co-chair.

Takeshi Naemura received his B.E., M.E., andPh.D. degrees in electronic engineering from the
University of Tokyo in 1992, 1994, and 1997. He
was a visiting associate professor of Stanford
University from 2000–2002. He is currently aprofessor of the Graduate School of Information
Science and Technology, the University of Tokyo.
His research interests include image-based rendering, 3D image processing, human interface,
mixed reality, and computational photography.