A Taxonomy and Evaluation of Dense Two-Frame Stereo Correspondence Algorithms

Daniel Scharstein of Middlebury College and Richard Szeliski of Microsoft Research published “A Taxonomy and Evaluation of Dense Two-Frame Stereo Correspondence Algorithms” in the International Journal of Computer Vision in 2002. A preliminary version, with Ramin Zabih as third author, appeared at the IEEE Workshop on Stereo and Multi-Baseline Vision (SMBV), held with CVPR 2001 in Kauai, and the full text was first released as Microsoft Research technical report MSR-TR-2001-81 in November 2001; this article follows the journal paper. The paper described dense stereo algorithms as combinations of four building blocks, measured those blocks one at a time on images with ground-truth disparities, and compared 20 complete algorithms. Its results table became the online Middlebury stereo evaluation. The task as understood today is covered in stereo matching, and the benchmark’s later history in Middlebury.

Problem

Stereo correspondence was one of the most heavily studied topics in vision, yet the authors found progress hard to gauge: most papers reported only qualitative results, and the last exhaustive surveys were about a decade old. In optical flow, the authors noted, Barron, Fleet, and Beauchemin’s 1994 comparison had had a dramatic effect: new methods now had to improve on published ones. Earlier stereo comparisons had evaluated sparse feature matchers, and Szeliski and Zabih’s 1999 workshop comparison, which this paper continued, had covered several dense algorithms without a taxonomy.

Comparing finished algorithms as black boxes, the authors argued, could not say whether a winning energy-minimization method owed its success to its energy or to its minimizer; answering that requires taking algorithms apart.

The scope was two rectified frames with known geometry and a disparity at every pixel, since view synthesis and image-based rendering need values even in occluded and textureless regions.

Contribution

• A taxonomy. Dense stereo algorithms perform subsets of four steps: matching cost computation, cost (support) aggregation, disparity computation or optimization, and disparity refinement. All operate on the disparity space image (DSI), the function over $(x, y, d)$ now usually called a cost volume. A table summarized 34 methods by matching cost, aggregation, and optimizer.
• A modular implementation. A stand-alone C++ test bed, released with the paper as StereoMatcher and no longer distributed, in which each step can be swapped and its parameters varied by script.
• Data with ground truth. New multi-frame image sets of piecewise planar scenes with subpixel-accurate disparities.
• An evaluation methodology. Error measures computed separately in regions where stereo typically fails.
• Experiments and a comparison. Controlled tests of each building block, and a comparison of 20 algorithms, 15 of them run by their own authors.

Method

The taxonomy

The four steps separate two families. Local (window-based) methods put their effort into the cost and its aggregation, then pick the lowest-cost disparity at each pixel (winner-take-all), with smoothness implicit in the window. Global methods usually skip aggregation and minimize an energy with an explicit smoothness term,

$$E(d) = \sum_{(x, y)} C(x, y, d(x, y)) + \lambda E_\text{smooth}(d),$$

and differ mainly in the minimizer: simulated annealing, mean-field diffusion, graph cuts, or dynamic programming on independent scanlines. Cooperative algorithms sit between the two, iterating local computations that behave like global optimization. Coarse-to-fine warping methods and multi-valued, voxel, and layered representations fell outside the framework.

The implementation

The test bed implemented basic versions of each component without extra post-processing. Costs were absolute or squared differences, optionally truncated and optionally made sampling-insensitive with Birchfield and Tomasi’s measure. Aggregation used box and binomial filters, shiftable windows (a box filter followed by a min-filter), and three diffusion schemes. Optimizers were winner-take-all, dynamic programming, simulated annealing, graph cuts using α-β swap moves, and a new scanline optimization (SO): dynamic programming along each row without ordering or occlusion constraints, solving the graph-cut problem minus its vertical smoothness terms. All global optimizers minimized the same energy, whose smoothness cost was multiplied by a penalty factor where the intensity gradient was small, so that disparity jumps would prefer intensity edges. Refinement fitted a parabola to the three costs around the winning disparity.

Data and measures

Each new sequence had 9 images from a camera on a translation stage. The scenes were built from posters and paintings, so hand-labeling planar regions and fitting each an affine motion gave ground-truth disparities. Of six sequences acquired, the study used two, Sawtooth and Venus. It added the University of Tsukuba “head and lamp” pair, whose hand-labeled ground truth has only integer disparities from 5 to 14, and the monochrome Map pair from Szeliski and Zabih’s study.

With computed disparities $d_C$ and ground truth $d_T$ over $N$ pixels, the evaluator computed the RMS error and the percentage of bad matching pixels,

$$R = \left( \frac{1}{N} \sum_{(x, y)} |d_C(x, y) - d_T(x, y)|^2 \right)^{1/2}, \qquad B = \frac{1}{N} \sum_{(x, y)} \big[ |d_C(x, y) - d_T(x, y)| > \delta_d \big],$$

with $\delta_d = 1$ pixel, chosen to coincide with earlier studies. Both were also computed over three regions derived from the reference image and ground truth, and their complements:

• Textureless pixels, where the squared horizontal intensity gradient, averaged over a small window, falls below a threshold.
• Occluded pixels, whose forward-mapped match lands where a nearer surface is visible.
• Depth discontinuity pixels, near jumps of more than 2 pixels in disparity, dilated by a $9 \times 9$ window.

Occluded pixels were excluded from the textureless and discontinuity statistics, and a border was ignored. The paper focused on three numbers, in its notation $B_{\bar{\mathcal{O}}}$ over non-occluded pixels, $B_{\bar{\mathcal{T}}}$ in textureless areas, and $B_{\mathcal{D}}$ near discontinuities, and preferred them to RMS error, which the few grossly wrong pixels dominate. It also defined a prediction error, comparing other frames of a sequence with images warped by the computed disparities.

Computing the Measures

The code below builds a synthetic rectified pair, a textured square at disparity 18 in front of a slanted background with a blank band, derives the three regions from the paper’s definitions and parameter values, and evaluates $5 \times 5$ and $15 \times 15$ SAD block matching with winner-take-all.

import numpy as np
from scipy import ndimage

rng = np.random.default_rng(0)
H, W, D = 120, 160, 24
y, x = np.mgrid[0:H, 0:W].astype(float)

# Ground truth (left image is the reference; left x matches right x - d):
# a slanted background plane and a fronto-parallel square in front of it.
gt = 4.0 + 0.04 * x
fg = (y >= 35) & (y < 85) & (x >= 70) & (x < 110)
gt[fg] = 18.0

# Texture lives on the surfaces; one band of the background is blank.
s = np.linspace(-40, W + 40, 4001)          # 20 samples per pixel
smooth = lambda n: ndimage.gaussian_filter1d(n, 20) / ndimage.gaussian_filter1d(n, 20).std()
tex_bg = 128 + 40 * smooth(rng.standard_normal(s.size))
tex_bg[(s > 15) & (s < 55)] = 128.0
tex_fg = 128 + 40 * smooth(rng.standard_normal(s.size))
row_bg = 20 * rng.standard_normal(H)        # vertical variation: an offset per row

def render(sees_fg, surf_x):
    """Sample the surface textures at surface coordinate surf_x and add sensor noise."""
    blank = ~sees_fg & (surf_x > 15) & (surf_x < 55)
    img = np.where(sees_fg, np.interp(surf_x, s, tex_fg),
                   np.interp(surf_x, s, tex_bg) + row_bg[:, None])
    return np.where(blank, 128.0, img) + rng.normal(0, 1.0, img.shape)

left = render(fg, x)
# Right image: pixel u sees the foreground at x = u + 18 if that is inside the square,
# otherwise the background point x with x - (4 + 0.04 x) = u.
xf = x + 18.0
sees_fg = (y >= 35) & (y < 85) & (xf >= 70) & (xf < 110)
right = render(sees_fg, np.where(sees_fg, xf, (x + 4.0) / 0.96))

def sad_wta(win):
    """Box-filtered absolute differences, winner-take-all, integer disparities."""
    cost = np.full((H, W, D), 255.0)
    for d in range(D):
        cost[:, d:, d] = np.abs(left[:, d:] - right[:, :W - d])
    cost = ndimage.uniform_filter(cost, size=(win, win, 1), mode="nearest")
    return cost.argmin(axis=2).astype(float)

# Evaluation regions, following the paper's definitions and Table 2 values.
gx2 = np.zeros((H, W))
gx2[:, 1:-1] = ((left[:, 2:] - left[:, :-2]) / 2) ** 2           # squared horizontal gradient
textureless = ndimage.uniform_filter(gx2, size=3) < 4.0           # eval_textureless width 3, thresh 4.0
# Occluded: some pixel further right lands left of this pixel's match in the right image.
u = x - gt
later_min = np.minimum.accumulate(u[:, ::-1], axis=1)[:, ::-1]
occluded = np.zeros((H, W), bool)
occluded[:, :-1] = later_min[:, 1:] < u[:, :-1]
jump = np.zeros((H, W), bool)                                     # eval_disp_gap = 2.0
jump[:, :-1] |= np.abs(np.diff(gt, axis=1)) > 2
jump[:-1, :] |= np.abs(np.diff(gt, axis=0)) > 2
discont = ndimage.binary_dilation(jump, np.ones((9, 9), bool))    # eval_discont_width = 9
valid = np.zeros((H, W), bool)
valid[10:-10, 10 + D:-10] = True                                  # border, plus full search range

regions = {
    "all": valid,
    "nonocc": valid & ~occluded,
    "occ": valid & occluded,
    "textured": valid & ~occluded & ~textureless,
    "textureless": valid & ~occluded & textureless,
    "discont": valid & ~occluded & discont,
}
print("pixels per region:", {k: int(m.sum()) for k, m in regions.items()})
for win in (5, 15):
    err = np.abs(sad_wta(win) - gt)
    print(f"\nSAD {win}x{win}, WTA   bad pixels B (>1 px)   RMS error R")
    for name, m in regions.items():
        print(f"  {name:12s} {np.mean(err[m] > 1.0):18.1%} {np.sqrt(np.mean(err[m] ** 2)):12.2f}")

Output:

pixels per region: {'all': 11600, 'nonocc': 11050, 'occ': 550, 'textured': 8984, 'textureless': 2066, 'discont': 1369}

SAD 5x5, WTA   bad pixels B (>1 px)   RMS error R
  all                       22.9%         4.88
  nonocc                    19.2%         4.26
  occ                       97.1%        11.70
  textured                   4.8%         2.19
  textureless               82.1%         8.73
  discont                   10.7%         3.34

SAD 15x15, WTA   bad pixels B (>1 px)   RMS error R
  all                       20.1%         5.02
  nonocc                    16.2%         4.29
  occ                       97.6%        12.71
  textured                   5.1%         2.67
  textureless               64.5%         8.21
  discont                   16.3%         4.04

The larger window cuts errors in the blank band from 82% to 65% but raises them near the square’s edges from 11% to 16%, the trade-off the paper found for every aggregation method. Occluded pixels are almost all wrong, which is why the paper kept them out of its main measure. In textured regions only about 5% of pixels are bad, yet the RMS error exceeds 2 pixels because those few errors are large.

Results

The component experiments, run mostly on Tsukuba, Sawtooth, and Venus, produced these findings:

• Matching costs. Absolute and squared differences performed about the same. Truncating the cost helped near discontinuities, typically at around 20 intensity levels, but truncating below the noise level made errors very large. With shiftable windows, untruncated costs did best, so the authors concluded that shifting the window to avoid outliers was preferable to limiting their influence.
• Aggregation. More aggregation helped in textureless regions and hurt near discontinuities, so the best extent was intermediate and image-dependent. Shiftable windows performed best, especially near discontinuities; the other methods performed similarly. The paper explained the “foreground fattening” of window methods: at a vertical depth edge, the object boundary itself is the strongest horizontal texture, and it matches at the foreground disparity.
• Optimization. Graph cuts consistently outperformed dynamic programming, scanline optimization, and simulated annealing, overall and in both problem regions, but took 10–30 minutes against under 2 seconds for the scanline methods. The best smoothness weight $\lambda$ varied by image, around 500 for Map and 20–200 for Tsukuba. Scanline methods showed horizontal streaks from their lack of inter-scanline consistency.
• Sub-pixel refinement. On a slanted newsprint region, parabola fitting lowered the RMS error of a $9 \times 9$ SSD matcher from 0.296 to 0.088 pixel, and matching at half-pixel steps before fitting brought it to 0.051. The authors recommended that combination to remove the staircase effect.

In the overall comparison, all methods ran with one parameter setting across the four image pairs and without sub-pixel estimation. Global methods on 2D Markov random fields generally performed best in all regions; most used graph cuts, and belief propagation also did well. Methods that modeled planar surfaces were especially strong on Sawtooth and Venus. Cooperative and diffusion methods often misplaced boundaries, and local and scanline methods did less well. On Tsukuba, the best $B_{\bar{\mathcal{O}}}$ was 1.15% (belief propagation), against 1.94% for the authors’ own graph-cut implementation. Per-image parameters improved some methods substantially, and the Map pair could be virtually “solved”. On speed, Hirschmüller’s real-time shiftable-window method took 0.1 seconds on Tsukuba with very good results; the authors’ graph-cut code took 662 seconds, against 23.6 seconds for Boykov and colleagues’ implementation with a newer max-flow algorithm.

Impact

The authors hoped a common dataset and methodology would become a “litmus test” for new stereo algorithms, and planned to maintain the comparison table online; 15 groups had sent disparity maps since the technical report’s release. Kendall and colleagues (2017) later credited the Middlebury and KITTI datasets with first motivating improved hand-engineered techniques for every component of stereo.

The test images saturated quickly: by 2003 Scharstein and Szeliski judged that progress was outpacing the ability of existing datasets to separate the best methods, and the Middlebury article traces the harder datasets that followed. Geiger, Lenz, and Urtasun (2012) introduced KITTI partly in reaction to the laboratory character of Middlebury, reporting that methods ranking high there performed below average on street scenes.

The taxonomy also outlived the methods it classified. Žbontar and LeCun’s MC-CNN presented its learned patch similarity as the matching-cost step of the four-step taxonomy and kept classical steps after it, and Kendall and colleagues’ GC-Net opened its related work with the same four steps before learning them end to end.

Limitations

• Simple scenes, little data. Piecewise planar scenes made precise ground truth possible but limited complexity, and the authors wrote that harder images would soon be needed. There were four pairs, one nearly solved, and Tsukuba’s ground truth was integer-valued over 10 disparities.
• Parameter sensitivity. Global methods depended strongly on $\lambda$, and dynamic programming on its occlusion cost, so fixed-parameter rankings partly reflect tuning.
• Basic components. The test bed deliberately omitted common post-processing such as left-right cross-checking and median filtering.

What Came After

Semi-global matching, which Hirschmüller likened to the paper’s scanline optimization run along many directions, is described in semi-global matching and in the paper article on its 2005 and 2008 publications. Stereo matching surveys the learned methods that followed, and cost volume the data structure they share.