Graph Cut based Continuous Stereo Matching using Locally Shared Labels
Tatsunori Taniai∗
University of Tokyo, Japan

Yasuyuki Matsushita
Microsoft Research Asia, China

Takeshi Naemura
University of Tokyo, Japan

taniai@iis.u-tokyo.ac.jp

yasumat@microsoft.com

naemura@nae-lab.org

Abstract

fields (MRFs) [8]. In terms of MRF optimization, however,BP is considered a sequential optimizer, which improves
each node individually keeping others conditioned at thecurrent state. In contrast, GC improves all nodes simultaneously by accounting for interactions across nodes, and thisglobal property helps optimization avoid local minima [20,25]. Nevertheless, incorporating spatial propagation into
GC-based optimization is not straightforward, because inference using GC proceeds rather all-nodes-simultaneously,
not one-by-one-sequentially like PatchMatch and BP.
In this paper, we introduce a new labeling scheme, locally shared labels, that enables spatial propagation infusion-based optimization using GC [17]. The locallyshared labels define, for each pixel or region, its compact
and local discrete label space that is shared among neighboring pixels/regions. By using locally shared labels wegenerate a number of disparity maps (so-called proposals in
the literature [17]), and fuse and refine them in an iterative
manner (see Fig. 1). For natural scenes that often exhibit locally planar structures, the joint use of locally shared labels
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 smooth solutions and to avoid
trapped at a bad local minima.
The advantages of our method are fourfold. First, ourlocally shared labels produce submodular moves that guarantee the optimal labeling at each min-cut (subproblem optimal), which in contrast is not guaranteed in general fusion
moves [17]. Second, this optimality property and spatial
propagation allow randomized search, rather than employ
external methods to generate plausible initial proposals as
done in previous fusion approaches [17, 25, 19], which may
limit the possible solutions. Third, our method achieves
greater accuracy than BP [3] thanks to the good properties
of GC and locally shared labels. Finally, unlike PMBP [3]the computation of both unary and pairwise costs can be
performed in a parallel manner1 , which is the most expensive part in practice. With the proposed approach, accurate

We present an accurate and efficient stereo matching
method using locally shared labels, a new labeling schemethat enables spatial propagation in MRF inference usinggraph cuts. They give each pixel and region a set of candidate disparity labels, which are randomly initialized, spatially propagated, and refined for continuous disparity estimation. We cast the selection and propagation of locallydefined disparity labels as fusion-based energy minimization. The joint use of graph cuts and locally shared labels
has advantages over previous approaches based on fusion
moves or belief propagation; it produces submodular moves
deriving a subproblem optimality; enables powerful randomized search; helps to find good smooth, locally planardisparity maps, which are reasonable for natural scenes;
allows parallel computation of both unary and pairwise
costs. Our method is evaluated using the Middlebury stereo
benchmark and achieves first place in sub-pixel accuracy.

1. Introduction
Recent years have seen significant progress in accuracy of stereo vision. One of the breakthroughs is the useof 3D labels [4, 3, 18, 19, 9]; by estimating a local 3D
disparity plane d = ap x + bp y + cp for each pixel, accurate photo-consistency is measured between matching pixels even with large matching windows. While stereo withstandard 1D discrete disparity labels [22, 15, 14, 6] can bedirectly solved by discrete optimizers such as graph cuts(GC) [16, 5] and belief propagation (BP) [26, 7], such approaches cannot be directly used for continuous 3D labels
due to the huge (infinite) label space (a, b, c) ∈ R3 .
Recent successful methods [4, 3, 18] use PatchMatch [1,2] to efficiently infer correct 3D planes using spatial propagation; each pixel’s candidate plane is, in raster-scan order,refined and then propagated to next pixels. Further in [3],this sequential algorithm is combined with BP yielding
an efficient optimizer PMBP for pairwise Markov random

1 Although BP is usually GPU-parallelizable, PMBP differs from BP’s
standard settings in that it defines label space uniquely and distinctively foreach pixel and propagate it; both make parallelization indeed non-trivial.

∗ Part of this work was done while the first author was visiting Microsoft

Research Asia as a research intern.

1

post-processing
random
1 iteration
2 iterations
3 iterations
10 iterations
Figure 1. Process of our stereo matching, showing disparity and 0.5-pixel error maps for Teddy. In our framework, we start with random
disparities that are represented by per-pixel 3D planes (leftmost). We then alternately propagate local planes using GC and refine them(middles). Finally, the resulting disparity map is further refined at post-processing stage based on left-right consistency check (rightmost).

stereo matching can be efficiently computed with a GPU
implementation as we will see in the experiment.

2. Related works
MRF stereo methods can be categorized into three approaches: discrete stereo, segment-based stereo, and continuous stereo.
Discrete stereo [22, 15, 14, 6] formulates stereo matching as a discrete multi-labeling problem, where each pixel
is individually assigned one of pre-defined discrete disparity values. For this problem, many powerful discrete optimizers, such as BP [26, 7], TRW [12], and GC [16, 5],
can be directly used. Successful results are shown using
GC with expansion moves [6, 20]. In expansion moves, themulti-labeling problem is reduced to a sequence of binarylabeling problems, each of which can be exactly solved by
GC, if only pairwise potentials ψ meet the following submodularity of expansion moves [13, 6]:ψ(α, α) + ψ(β, γ) ≤ ψ(β, α) + ψ(α, γ).

(1)Segment-based stereo [21, 10, 11, 23] assigns a 3D disparity plane for each of over-segmented image regions.
The candidate planes are generated by fitting planes to aroughly estimated disparity map, and then the optimal assignment of the planes is estimated by, e.g., GC with expansion moves [6, 10] or BP [7, 11]. Although this approachyields continuous-valued disparities, it strictly limits the reconstruction to a piecewise planar representation. Also, results are subject to the quality of the segmentation.
The last group, to which our method belongs, is continuous stereo [25, 4, 3, 19, 18, 9], where each pixel is assigned a distinct continuous disparity value. Some methods [25, 19] use fusion moves [17], an operation that combines two disparity maps to make a better one (binary fusion) by solving a non-submodular binary-labeling problem using QPBO-GC [13, 17]. In this approach, a numberof continuous-valued disparity maps (or proposals) are first

generated by other external methods (e.g., segment-basedstereo [25]), which are then combined as a sequence of binary fusions. Our method is also based on fusion moves
but generates proposals using locally shared labels, whichenable spatial propagations of local candidate planes and,more importantly, they make fusion moves submodular, i.e.,each binary fusion is optimally solved via GC (subproblem
optimal). Our method only requires randomized initial proposals instead of those generated by external methods. A
stereo method by Bleyer et al. [4] proposes accurate photoconsistency measures using 3D disparity planes that are inferred by PatchMatch [1, 2]. Heise et al. [9] incorporates
Huber regularization into [4] using convex optimization.
Besse et al. [3] point out a close relationship between PatchMatch and BP and present a unified method called PatchMatch BP (PMBP) for pairwise continuous MRFs. PMBP
is probably the closest approach to ours 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 [20]for achieving greater accuracy. In addition, our method allows parallel computation of both unary and pairwise costs.

3. Proposed method
This section describes the proposed stereo matching
method. Given two input images IL and IR , our purpose
is to estimate disparity maps of both images.

3.1. Formulation
We use a pairwise MRF formulation by following conventional stereo matching methods [19, 22, 15, 14, 6]. In
the MRF framework, each pixel p ∈ (P ⊂ Z2 ) is assigned
a value in some disparity space S, and one seeks a disparity
map f for every pixel fp = f (p) : P → S that minimizes
XX X
E(f ) =
φp (fp ) + λψpq (fp , fq ).
(2)
p∈Pp∈P q∈N (p)The first term, called the data term or unary term, measures the photo-consistency between matching pixels. The

disparity 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 of a pixel p and its neighboring pixels q ∈ N (p). We define these terms as below.
Data term. To measure photo-consistencies, we use adata term that has been recently proposed by [4]. Here,each pixel p’s disparity dp is over-parameterized by a 3D
plane dp = ap x + bp y + cp to avoid the frontal-parallelbias. Therefore, the objective becomes to seek a disparity plane fp = (ap , bp , cp )T ∈ S for every pixel p in the left
and right images such that disparity map f minimizes theenergy function E(f ) of Eq. (2). Using this p’s disparity
plane fp , a pixel q = (qx , qy )T in the left image is warped
to a new location in the right image by a warping wfp as
wfp (q) = q − (ap qx + bp qy + cp , 0)T .

(3)The data term of p in the left image is therefore defined asXφp (fp ) =ωpq ρ q, wfp (q) .
(4)
q∈Wp

Here, Wp is a square window centered at p. The weight ωpq
implements the adaptive support window proposed in [27],and is defined asωpq = e−kIL (p)−IL (q)k1 /γ ,(5)where γ is a user-defined parameter, and k · k1 represents
the `1 -norm. The function ρ(q, wfp (q)) measures the pixel
dissimilarity between q and its matching point wfp (q) asρ(q, wfp (q)) = (1−α) min(kIL (q)−IR (wfp (q))k1 , τcol )+α min(k∇x IL (q)−∇x IR (wfp (q))k1 , τgrad ), (6)where ∇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 to compute
IR (wfp (q)). When the data term is defined on the right image, we swap IL and IR in Eqs. (5) and (6), and add thedisparity value in Eq. (3).
Smoothness term. For the smoothness term, we use
a curvature-based, second-order smooth regularization
term [19] defined as
ψpq (fp , fq ) = max (ωpq , ) min(ψ̄pq (fp , fq ), τdis ), (7)where  is a small constant value that gives a lower bound to
the weight ωpq for increasing the robustness. 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 )|,(8)

where dp (fq ) = aq px + bq py + cq . ψ̄pq (fp , fq ) is truncated by τdis to allow sharp jumps in disparity at depthedges. The term ψpq (fp , fq ) satisfies the submodularity ofEq. (1) for taking advantage of GC. Visualization and proof
for submodularity of this term are shown supplementary.

3.2. Locally shared labels
As the main contribution of this paper, we introduce
locally shared labels for efficiently optimizing continuous
MRFs. The locally shared labels are the combination of
pixel and region labels, in which label spaces are shared
among neighbors, and they enable per-pixel estimation ofcontinuous solutions as well as fast propagations.
Pixel and region labels. Pixel labels are a small number(say, K) of discrete disparity labels (or candidate labels) defined at each pixel p, which we refer to as a pixel label set,(0) (1)(K−1)(i)
Lp = {lp , lp , . . . , lp}, lp = (a, b, c)T ∈ S. The
pixel label sets are shared among neighboring pixels. In
addition, we define region labels that give additional candidate labels for accelerating spatial propagation and avoiding stuck at a local minima. We use a regular grid structure
for regions, which are indexed by the region coordinatesr ∈ (R ⊂ Z2 ) like the pixel coordinates. Region labels
define for a region r a set of KR candidate labels Rr ⊂ S,
which we call a region label set. Each label set Rr gives
candidate labels for pixels in the region r and pixels in theneighboring regions as well, i.e., Rr is also shared amongneighboring regions just like pixel labels.
During the inference, for each pixel p, our methodchooses the best candidate label fp from the union of pixel
and region label sets that are shared for the pixel p:!
!
[[
Cp = Lp ∪Lq ∪ Rr ∪
Rs ,
(9)
q

s

where q, r, s represent p’s neighboring pixels, the region
that p belongs to, and the neighboring regions to r, respectively. By sharing local label sets among neighbors, goodcandidate labels are spatially propagated to nearby pixels.
The concept of pixel and region labels is illustrated in Fig. 2.
Proposal generation for fusion. During the inference,
we repeatedly seek the best labeling f (t) for the current local label sets {Lp } and {Rr }, and refine them. The former part, i.e., the selection and propagation of candidate
labels in {Lp } and {Rr }, is cast as fusion-based energy
minimization [17], as described in the rest of this section.
Consider the essential function of fusion is to make a goodsolution f by fusing a number of proposal disparity maps
{g (0) , g (1) , . . . , g (n−1) }, where at each pixel fp is assigned(0) (1)(n−1)
one of n disparity labels {gp , gp , . . . , gp}. In our
method, we build a special form of proposals from {Lp }and {Rr } in a manner that achieves the propagation of pixel
and region labels. To make proposals from pixel labels, wetop view of
(0)

�
�(1)

�0

�5

�2

�1

�3

(K = 1)

�2

�6(a) Pixel labels (1D image)sharedregion�3�7(b) Region labels (1D image)

…

�(2)
�(3)

�0

�4
�1

proposals(c) Proposals in case of pixel labels

Figure 2. Illustrations of pixel and region labels, and proposal construction. For simplicity, they are illustrated by 1D images with the blue
nodes representing pixels. The boxes Lp and Rr represent a set of candidate disparity labels given for the pixel p and pixels in the regionr, respectively. The horizontal red lines signify that the local label sets {Lp } and {Rr } are shared between neighbors. The local label sets
are aligned so as to make proposal disparity maps g (j) for fusion. In particular, proposal construction of pixel labels is illustrated for a 2D
image in (c) indicating that, in each proposal g (j) , candidate labels defined at blue pixels are shared between red neighbors.
(i)

copy each candidate label lp ∈ Lp to the j-th proposal(j)(i)
g (j) by setting gq ← lp , where q is the nine neighboring pixels around and including p (i.e., the pixels where the(i)candidate label lp is shared), and j is given as
j = K (4 (py mod 4) + (px mod 4)) + i.

(10)Here, mod is the modulo operation, and px and py
are p’s coordinates specified as px ∈ [0, width − 1] andpy ∈ [0, height − 1]. Figure 2a illustrates this construction(0)
for the case of a 1D image with K = 1 (i.e., Lp = {lp })for simplicity, where a horizontal layer of candidate labels
at vertical position j represents a proposal g (j) . Figure 2cillustrates for a 2D image showing that candidate labels atblue pixels are shared among red neighbors in each proposal g (j) . The integer 4 in Eq. (10) means that, in each
proposal, we leave a “gap” (shown as gray pixels in Fig. 2cthat represent no candidate labels) between each “shared region” (see Fig. 2c) for ensuring submodularity, which we
describe later. We assign an infinite unary cost to those invalid labels to ensure that such labels are avoided during theinference. For region labels, proposals are constructed in
the same manner with pixel labels by regarding a region asa pixel as shown in Fig. 2b. The fusion is performed using
the proposals generated from both pixel and region labels.
This particular form of proposal construction guarantees
that a binary fusion of an arbitrary solution f and any ofthe proposals g = g (j) is submodular, thus it is exactly
solved via GC. To show this, we give an intuitive explanation here due to the limited space. A complete proof is insupplementary. Let’s consider a simple case where gp takes
the same value for all pixels p. A fusion move with suchglobally-constant g is called an expansion move [6], which
is submodular if pairwise terms ψpq satisfy Eq. (1). Seeing
Fig. 2c, our proposal g is made locally-constant by “sharedregions”, thereby a fusion move with our proposals is, virtually, many mutually-disjoint local-alpha-expansions withdifferent “alpha” for each shared region. Therefore the condition for deriving submodularity is the same with Eq. (1),which holds for our smoothness term ψpq of Eq. (7). With

Algorithm 1 Overview of optimization procedure
1: Initialize {Lp } and {Rr } randomly.
2: repeat3:♦ Optimize labeling f for current local label sets:

f (t) = argmin E(f ) with local label sets {Lp } and {Rr }4:♦ Refine local label sets {Lp } and {Rr }:5:
for all pixels p ∈ P do
6:C̃p ← Cp with perturbation.
(t)
7:
Lp ← best K−1 candidate labels c ∈ (Lp ∪C̃p )\{fp }(t)that minimize Ep (c|f )(t)
8:
Lp ← Lp ∪ {fp }9:end for10:
for all regions r ∈ R do(t)
11:
Rr ← random KR candidate labels from {fp |p in r}
12:
end for13: until convergence

this submodularity guarantee, we only need to use a standard GC [16, 5] instead of employing expensive QPBOGC [13] used in usual fusion approaches [17, 19, 25].
In addition, this proposal generation helps obtain smooth
solutions because multiple pixels in shared regions are allowed to move-at-once to the same candidate label at one
binary fusion. This effect becomes more significant with
region labels because of their large shared regions. In fact,region labels make the key factor in our algorithm for both
efficiency and accuracy as we will see in the experiment.

3.3. Optimization
The overview of our optimization procedure is summarized in Algorithm 1. As discussed in the previous section,
our optimization uses an iterative framework, where we alternately optimize the labeling f with given local label sets{Lp } and {Rr }, and refine the local label sets {Lp } and{Rr } locally with the labeling f fixed.
Our optimization begins with randomly initializing {Lp }and {Rr }. To sample the allowed solution space evenly,
we take the initialization strategy described in [4]. For(i)
lp ∈ Lp at p = (px , py )T , we select a random disparity z0

in the allowed disparity range [0, dispmax]. Then, a random unit vector n = (nx , ny , nz )T and z0 are converted tothe plane representation by ap = −nx /nz , bp = −ny /nz ,and cp = −(nx px + ny py + nz z0 )/nz . For the region label sets {Rr }, we randomly pick KR pixels in each region,(0)
and copy the candidate label lp ∈ Lp of the randomly chosen pixels p to the region label set.
At line 3 of Algorithm 1, we optimize the labeling f by
the procedure illustrated in the previous section. It can beapproximately solved by sequentially fusing proposals constructed from pixel and region label sets {Lp } and {Rr }.
In lines 5–9, we refine the pixel label sets {Lp }. At eachpixel p, we first randomly perturb p’s candidate labels Cp of
Eq. (9) and obtain C̃p . As the refined Lp , we select the best
K candidate labels from the union of C̃p and the current Lpthat minimize the following local energy at the pixel p:Xψpq (s, fq(t) ) , s ∈ S. (11)Ep (s|f (t) ) = φp (s) +q∈N (p)Here, the refined Lp is forced to contain the current candi(t)date label fp to ensure that, in the next iteration, the solution f (t+1) can stay at f (t) , thereby the energy does notincrease, i.e., E(f (t) ) ≥ E(f (t+1) ) holds throughout the iterations. Perturbation is implemented as described in [4].
Namely, each candidate label (a, b, c)T ∈ Cp is converted
to disparity d and normal vector n. We then add a random
disparity ∆d ∈ [−rd , rd ] and a random unit vector ∆n tothem, respectively, as d0 = dp + ∆d and n0 = n + rn ∆n .
Finally, d0 and n0 /|n0 | are converted to the plane representation (a0 , b0 , c0 )T ∈ C̃p as a perturbed candidate label. The
values rd and rn define an allowed change of planes. We
start by setting rd ← dispmax/2 and rn ← 1. After each
iteration, we update them by rd ← rd /2 and rn ← rn /2.
In lines 10–12, we update the region label sets {Rr }. As
done in the initialization, we again take a random-pick-up(t)scheme. This time, the current solution fp of randomly
chosen pixels p is taken as the region labels.
Finally, after the whole process, we perform the postprocessing using left-right consistency check and medianfiltering as described in [4] for further improving the results.
This step is widely employed in recent methods [4, 3, 18, 9].

4. Experiments
In the experiments, we first evaluate our method on the
Middlebury benchmark. We further assess the effect of region labels, and also compare with the PMBP method [3]that is closely related to our approach.
Setting. We use the following settings throughout the experiments. We use a PC with a Xeon CPU (2.53 GHz ×
4 cores) and NVIDIA GeForce GTX-295 GPU (only oneof two GPU cores is used here). The parameters of ourdata term are set as {τcol , τgrad , γ, α} = {10, 2, 10, 0.9}

as specified in [4]. The size of supporting windows is set
to 41 × 41, which is the same setting as PMBP [3]. For thesmoothness term, we use {λ, τdis , } = {20, 1, 0.01} andeight neighbors for N . For optimization, three-layer locallyshared labels are used: pixel labels with K = 2, region labels of size 5 × 5 with KR = 2, and also regions labels of
size 25 × 25 with KR = 2, and a GC implementation of [5]is used. We iterate twice for each proposal in fusion stage,and iterate the outer-loop process ten times. The computation of unary costs is performed in parallel on GPU, andpairwise costs are computed on four CPU cores.

4.1. Evaluation on the Middlebury benchmark
We show in Tab. 1 selected rankings on the Middlebury stereo benchmark for 0.5-pixel accuracy. Our method
achieves the current best average rank (3.5) and bad-pixelrate (6.63%) amongst more than 155 stereo methods. Even
without post-processing, our method still outperforms theother methods in average rank, despite that methods [4, 3,
18, 9] use the post-processing. Compared with closely related approaches (PMBP [3] and PatchMatch stereo [4]),which are ranked seventh and ninth in Tab. 1, although
results of PMBP for Cones are slightly better than ours,
our method consistently outperforms the two methods inthe other evaluations. We summarize the results of our
method in Fig. 3. We omit the results of Tsukuba, because
its ground truth is quantized to integer disparities; thus weconsider that not to be appropriate for sub-pixel evaluation.
More results including failures are shown in supplementary.

4.2. Effect of region labels
To observe the effect of region labels, we assess the performance using three different settings: (1) only pixel labelswith K = 6; (2) pixel labels with K = 4 and region labelsof size 5 × 5 with KR = 2; (3) pixel labels with K = 2,
region labels of size 5 × 5 with KR = 2, and also regions
labels of size 25 × 25 with KR = 2 (the default setting forour method described above). Among the three settings, thenumber of candidate labels given for each pixel (i.e., |Cp |)is kept consistent. We use λ = 40 but keep the other parameters as default. Using these settings, we observe the performance variations by estimating the disparities of only theleft image of the Cloth1 dataset without the post-processing.
In Fig. 4, we show the results of the three cases after ten
iterations. Also, plots in Figs. 5a-c show the energy variations of the energy function, data term, and smoothness term
after each iteration, respectively. Figure 5a shows that region labels play a critical role in minimizing energies. Figures 5b and 5c indicate that region labels effectively reduce
the energies of both the data and smoothness terms, whichshows the contribution of region labels for fast propagationof good candidate labels as we intended. On the other hand,
if only pixel labels are used, the solution is trapped at a badTable 1. Middlebury benchmark for 0.5-pixel accuracy. Our method achieves the current best average rank 3.5. In all, results are evaluatedfor all pixels where the ground truth is given, while only for non-occluded pixels in nonocc, and around depth discontinuities in disc.

Image
After post-processing
0.5-pixel error map
Raw result
Ground truth
Figure 3. Our results on the Middlebury benchmark. From left to right, one of the input images, our results without post-processing, after
post-processing, the ground truth, and 0.5-pixel error maps of the results after post-processing are shown. In the error maps, white andblack pixels indicate correct and incorrect disparities, while gray indicates incorrect but occluded pixels.

local minima producing a noisy result as in Fig. 4b. Although the use of the large region labels does not make asignificant difference in the converged energy values, thedifference is obvious if we see Figs. 4c and 4d.

4.3. Comparison with PMBP [3]We compare our method with PMBP [3] that is the closest method to ours. For a fair comparison, we use fourneighbors for N in Eq. (2), which is the same setting as
PMBP. For a comparable smoothness weight with the default setting (eight-neighbor N ), we use λ = 40 and keep
the other parameters as default. For PMBP, we use the samemodel as ours; the only difference from the original PMBP
is the smoothness term, which does not satisfy the submodularity of Eq. (1). PMBP also defines K candidate labels
for each pixel, for which we set K = 1 and K = 5 (origi-

nal paper uses K = 5). We show the comparison using theCones dataset by estimating the disparity map of only theleft image without the post-processing.
Figures 6a-c show the temporal transition of the energy
values in the full and zoomed scales, and the 0.5-pixel error
rates, respectively. We show the performance of our method
using its GPU and CPU (1 or 4 CPU-cores) implementations. For PMBP, we also implemented the unary cost computation on GPU, but it became rather slow, possibly due tothe overhead of data transfer. Efficient GPU implementations for PMBP are not available in literature2 . Therefore,
2 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 it propagatescandidate labels to distant pixels, it is not applicable to PMBP that must
propagate messages to neighbors, and is not as efficient as our 100x, either.

(c) Pixel and region(d) Pixel and region(e) Ground truth(5 × 5) labels(5 × 5 and 25 × 25) labels
Figure 4. Visual effect of region labels. (a) One of input images. (b) Using only pixel labels yields a noisy result, which is improved by (c)adding region labels. (d) Large region labels are effective for occluded regions. These are all raw results without the post-processing.

Pixel labels

Pixel + 5x5 region labels

Pixel +
5x5 region +25x25 region labels
0

2

4
6
# iterations

8

10

3.7
3.6
3.5
3.4
3.3
3.2
3.1
3.0
2.9
2.8
2.7
2.6
2.5

14.0
Smoothness term value [105 ]Energy function value [107 ]

3.7
3.6
3.5
3.4
3.3
3.2
3.1
3.0
2.9
2.8
2.7
2.6
2.5

Data term value [107 ](b) Only pixel labels(a) Cloth1

Pixel labels
Pixel + 5x5 region labels
Pixel +
5x5 region +25x25 region labels
02

Pixel labels

12.0

Pixel + 5x5 region labels

10.0
8.0

Pixel + 5x5 region +25x25 region labels

6.0
4.0
2.0

0.0

4
6
# iterations

8

10

0

2

4
6
# iterations

810(a) Energy function(b) Data term(c) Smoothness term
Figure 5. Effect of region labels in minimizing energies. From left to right, energy transitions of the energy function, data term, andsmoothness term w.r.t. the number of iterations are shown. Region labels play a critical role in reducing the energies.

the plots show PMBP results that use a single CPU core.
Figures 6a and 6b show that PMBP works much faster thanour CPU implementation; however, our GPU implementation shows significantly faster convergence. Furthermore,
our method reaches the better solution than that of PMBPin both energy values and error rates3 . At around 4000[sec]of Figs. 6b and 6c, the solution obtained by our CPU implementation marked the lower error rate than that of PMBPin spite of its higher energy. Figure 7 shows the resulting disparity maps obtained by our method and PMBP with
K = 5. Our result shows greater accuracy around the edge
regions.

5. Conclusions
In this paper, we presented an accurate and efficientstereo matching method for continuous disparity estimation. Unlike previous approaches that use fusion [17,19, 25], our method is subproblem optimal and only requires randomized initial proposals. By comparing witha recent continuous MRF stereo method, PMBP [3], ourmethod showed an advantage in efficiency and comparable
or greater accuracy. The use of a GC-based optimizer makes
our method advantageous.
3 The CPU implementation of our method also reaches almost the sameenergy and error rate after about 42000[sec] by a single core, 9150[sec]by four cores.

We believe that our optimization method can be applied
for more general corresponding field estimation such as optical flow, but we leave it for our future work. We also
believe that some occlusion handling schemes based on
GC optimization [14, 15, 24] can be incorporated into ourframework, which may yield even greater accuracy.

References[1] 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. 1, 2, 7
[2] C. Barnes, E. Shechtman, D. B. Goldman, and A. Finkelstein. The
generalized patchmatch correspondence algorithm. In Proc. of European Conf. on Computer Vision (ECCV), pages 29–43, 2010. 1,
2[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. 1,
2, 5, 6, 7, 8
[4] 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. 1, 2, 3, 5
[5] Y. Boykov and V. Kolmogorov. An experimental comparison of mincut/max- flow algorithms for energy minimization in vision. IEEE
Trans. Pattern Anal. Mach. Intell., 26(9):1124–1137, 2004. 1, 2, 4, 5
[6] 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. 1, 2, 4

3.2

PMBP
(K=1)

Ours
(CPUx1)
PMBP
(K=5)

Ours
(GPU)

26

PMBP
(K=1)

3.1

PMBP
(K=5)

Ours
(CPUx4)

3.0
Ours
(GPU)

2.9
0

1000

2000 3000 4000
Time [seconds]

5000

Sub-pixel error rate [%]Ours(CPUx4)

Energy function value [107 ]Energy function value [107 ]

4.2
4.1
4.0
3.9
3.8
3.7
3.6
3.5
3.4
3.3
3.2
3.1
3.0
2.9

24
22
Ours
(CPUx4)

PMBP
(K=1)

20
18
16
14

Ours
(GPU)

12
10

PMBP
(K=5)

8
0

1000

2000 3000 4000
Time [seconds]

5000

0

1000

2000
3000
4000
Time [seconds]

5000(a) Time-energy transition (full)(b) Time-energy transition (zoomed)(c) Error rate
Figure 6. Efficiency and accuracy comparison with PMBP [3]. Accuracies are evaluated for all-regions after each iteration.

(a) Our method(b) PMBP [3](c) Ground truth
Figure 7. Visual comparison with PMBP. Our method finds better disparities around edges and occluded regions without post-processing.

[7] P. Felzenszwalb and D. Huttenlocher. Efficient belief propagation forearly vision. In Proc. of IEEE Conf. on Computer Vision and Pattern
Recognition (CVPR), volume 1, pages 261–268, 2004. 1, 2
[8] 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. 1[9] 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. 1, 2, 5
[10] L. Hong and G. Chen. Segment-based stereo matching using graphcuts. In Proc. of IEEE Conf. on Computer Vision and Pattern Recognition (CVPR), volume 1, pages 74–81, 2004. 2[11] A. Klaus, M. Sormann, and K. Karner. Segment-Based Stereo
Matching Using Belief Propagation and a Self-Adapting Dissimilarity Measure. In Proc. of Int’l Conf. on Pattern Recognition (ICPR),
volume 3, pages 15–18, 2006. 2[12] V. Kolmogorov. Convergent Tree-Reweighted Message Passing forEnergy Minimization. IEEE Trans. Pattern Anal. Mach. Intell.,28(10):1568–1583, 2006. 2[13] V. Kolmogorov and C. Rother. Minimizing Nonsubmodular Functions with Graph Cuts – A Review. IEEE Trans. Pattern Anal. Mach.
Intell., 29(7):1274–1279, 2007. 2, 4
[14] V. Kolmogorov and R. Zabih. Computing visual correspondence withocclusions using graph cuts. In Proc. of Int’l Conf. on Computer
Vision (ICCV), volume 2, pages 508–515, 2001. 1, 2, 7
[15] V. Kolmogorov and R. Zabih. Multi-camera Scene Reconstruction
via Graph Cuts. In Proc. of European Conf. on Computer Vision(ECCV), pages 82–96, 2002. 1, 2, 7
[16] 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. 1, 2, 4
[17] V. Lempitsky, C. Rother, S. Roth, and A. Blake. Fusion Movesfor Markov Random Field Optimization. IEEE Trans. Pattern Anal.
Mach. Intell., 32(8):1392–1405, 2010. 1, 2, 3, 4, 7

[18] 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. 1, 2, 5
[19] 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. 1, 2, 3, 4, 7
[20] 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. 1, 2
[21] 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. 2[22] 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. 1, 2
[23] Z.-F. Wang and Z.-G. Zheng. A region based stereo matching algorithm using cooperative optimization. In Proc. of IEEE Conf. on
Computer Vision and Pattern Recognition (CVPR), 2008. 2[24] 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. 7[25] 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. 1, 2, 4, 7
[26] 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. 1, 2
[27] K.-J. Yoon and I.-S. Kweon. Locally adaptive support-weight approach for visual correspondence search. In Proc. of IEEE Conf. on
Computer Vision and Pattern Recognition (CVPR), volume 2, pages
924–931, 2005. 3