Stereo Rectification & Disparity Estimation

"I compare a million pixels against a million pixels before breakfast, and all I ask in return is that somebody straightens the epipolar lines first. I match along rows. I do not do diagonals. Have some respect."

An Uncompromisingly Horizontal Block Matcher

Everything so far in this chapter operated on sparse keypoints: a few hundred matches, carefully verified. Applications want more. A robot avoiding obstacles, a car braking for a pedestrian, a phone faking background blur: all need depth at every pixel, which means correspondence at every pixel, including the boring, textureless ones that Chapter 10's detectors deliberately skip. Dense matching is only computationally sane because of the epipolar constraint, and it becomes downright pleasant when the constraint is given its friendliest shape.

1. Rectification: Engineering the Easy Geometry Intermediate

Section 13.1 identified the dream configuration: cameras side by side, optical axes parallel, epipoles at infinity, epipolar lines horizontal and row-aligned. Physical rigs only approximate it; cameras are often slightly verged (angled inward toward each other so their fields of view overlap sooner), and lens distortion and millimeter-level mounting errors leave corresponding points a few rows apart, which is fatal for a matcher that searches only along rows. Rectification fixes this in software: rotate each camera (virtually, by warping its image) about its own optical center until the dream configuration holds exactly. Rotating about the center changes no ray origins, so this is the rotation-only homography of Section 13.3 applied to each image, computed so that both new image planes are coplanar and parallel to the baseline. Figure 13.4.1 shows the before and after.

With calibration from Chapter 12 (intrinsics of both cameras plus the rig's rotation and translation R, T from stereo calibration), OpenCV computes the rectifying rotations and the new projection matrices in one call, and the warps in two more. These maps are computed once per rig and reused for every frame:

# One-time rectification setup for a calibrated stereo rig: turn the rig's
# intrinsics and relative pose into per-pixel warp maps, then verify the
# warp by checking that matched features land on nearly the same row.
import cv2
import numpy as np

# K1, dist1, K2, dist2, R, T come from cv2.stereoCalibrate (Chapter 12)
R1, R2, P1, P2, Q, roi1, roi2 = cv2.stereoRectify(
    K1, dist1, K2, dist2, (w, h), R, T,
    alpha=0)                      # alpha=0: crop to valid pixels only

mapLx, mapLy = cv2.initUndistortRectifyMap(K1, dist1, R1, P1, (w, h), cv2.CV_32FC1)
mapRx, mapRy = cv2.initUndistortRectifyMap(K2, dist2, R2, P2, (w, h), cv2.CV_32FC1)

rectL = cv2.remap(imgL, mapLx, mapLy, cv2.INTER_LINEAR)   # per-frame, fast
rectR = cv2.remap(imgR, mapRx, mapRy, cv2.INTER_LINEAR)

# sanity check: matched features should now land on (nearly) equal rows;
# on a good calibration the mean row difference is below 0.5 pixels
sift = cv2.SIFT_create(nfeatures=500)
kL, dL = sift.detectAndCompute(rectL, None)
kR, dR = sift.detectAndCompute(rectR, None)
mm = [m for m, n in cv2.BFMatcher(cv2.NORM_L2).knnMatch(dL, dR, k=2)
      if m.distance < 0.7 * n.distance]
dy = [abs(kL[m.queryIdx].pt[1] - kR[m.trainIdx].pt[1]) for m in mm]
print(f"mean |row difference|: {np.mean(dy):.2f} px")   # e.g. 0.31 px

Two practical notes. The alpha parameter trades field of view against black borders: 0 crops to fully valid pixels, 1 keeps everything. And the returned $Q$ matrix is not rectification's business at all; it is the disparity-to-depth reprojection that Section 13.5 runs on, so keep it. For pairs without calibration, cv2.stereoRectifyUncalibrated computes rectifying homographies from a fundamental matrix (Section 13.2), good enough for visualization though not for metric depth.

2. Disparity, and the Simplest Matcher That Could Work Beginner

On a rectified pair, a scene point projects to $(x_L, y)$ on the left and $(x_R, y)$ on the right, and since the right camera sits to the right, $x_R \le x_L$: everything shifts left, near things more than far things (hold a finger close to your face and alternate eyes). The disparity $d = x_L - x_R \ge 0$ is the entire output of stereo matching: a map of horizontal shifts, large for near pixels, small for far ones, inversely proportional to depth in a way Section 13.5 makes exact.

The simplest dense matcher is block matching: for each left pixel, slide a small window along the same row of the right image over a disparity range $[0, d_{\max}]$, score each shift by a dissimilarity such as the sum of absolute differences (SAD), and keep the minimizer (the illustration below shows the slide in action):

$$ d(x, y) \;=\; \arg\min_{d' \in [0, d_{\max}]} \sum_{(u,v) \in W} \big| I_L(x+u,\, y+v) - I_R(x+u-d',\, y+v) \big|. $$

A single pixel carries too little information to match: an 8-bit pixel takes one of 256 gray values, so on a 1000-pixel row roughly four pixels share each value, and matching a lone pixel means choosing among them by coin flip. Hence the window $W$: a $7 \times 7$ patch of those same pixels has $256^{49}$ possible appearances, astronomically more than there are pixels to confuse it with, so the matching unit becomes a neighborhood, on the bet that neighborhoods are distinctive where pixels are not. That bet, and where it loses, structures the rest of the section.

Block matching fails in three canonical places: textureless regions (blank walls: every shift scores alike, the minimum is noise), repetitive structure (fences, brick, racking: several shifts score well, the wrong minimum wins), and occlusions (pixels near depth edges visible to only one camera: the true match does not exist, but the matcher returns a minimum anyway). A fourth, subtler failure is edge fattening: a window straddling a depth boundary contains two motions, and the textured side wins, bleeding foreground disparity into the background.

3. Semi-Global Matching: Smoothness at Scanline Prices Advanced

Semi-global matching (SGM, Hirschmüller 2008) keeps a per-pixel matching cost $C(\mathbf{p}, d)$ but adds a smoothness penalty: a small constant $P_1$ for neighbors whose disparity differs by exactly 1 (sloped surfaces are fine), a larger $P_2$ for jumps greater than 1 (depth discontinuities are allowed but must pay). Minimizing the joint energy over the full 2D grid is NP-hard, but along a single 1D path it is classic dynamic programming. SGM's move: run the 1D optimization along 8 (or 16) directions through each pixel (horizontal, vertical, diagonals), sum the directional costs, and pick each pixel's best disparity from the aggregate. Smoothness information flows in from all around the image at a cost linear in pixels times disparities, and the streaking artifacts of single-direction scanline stereo cancel out across paths.

OpenCV's StereoSGBM ("semi-global block matching") implements a practical variant whose handful of parameters do real work, annotated below. The non-obvious convention: compute returns disparities as 16-bit integers scaled by 16, so divide before use.

# Configure semi-global block matching: each parameter encodes a piece of
# the SGM theory (the P1/P2 smoothness penalties, the uniqueness and
# left-right checks), and compute() returns disparity as int16 scaled by 16.
block = 5
stereo = cv2.StereoSGBM_create(
    minDisparity=0,
    numDisparities=128,            # search range; MUST be divisible by 16.
                                   # ceiling on measurable nearness (Sec 13.5)
    blockSize=block,               # odd; 3-9: small=detail, large=smooth
    P1=8 * 3 * block ** 2,         # penalty for |delta d| == 1 (slopes)
    P2=32 * 3 * block ** 2,        # penalty for jumps; P2 > P1 enforces
                                   # piecewise-smooth surfaces
    uniquenessRatio=10,            # best cost must beat runner-up by 10%
    speckleWindowSize=100,         # remove disparity blobs smaller than this
    speckleRange=2,                # ...if they jump by more than 2 levels
    disp12MaxDiff=1,               # left-right consistency tolerance (px)
    mode=cv2.STEREO_SGBM_MODE_SGBM_3WAY)   # faster 3-path aggregation

disp = stereo.compute(rectL, rectR).astype(np.float32) / 16.0
valid = disp > 0                   # negative/zero = no reliable match
print(f"valid: {100 * valid.mean():.1f}%, "
      f"range: {disp[valid].min():.1f} to {disp[valid].max():.1f} px")
# valid: 91.3%, range: 4.2 to 97.8 px

The quality-control parameters deserve emphasis because they implement the failure analysis of subsection 2. uniquenessRatio attacks repetitive structure: if the runner-up disparity scores nearly as well as the winner, the pixel is marked invalid instead of guessed. disp12MaxDiff enforces left-right consistency: matching right-to-left as well, and invalidating pixels whose two answers disagree, which catches most occluded pixels (their two matches land on different surfaces). The speckle filter removes small floating islands of disparity, usually mismatches in low texture. A mature pipeline often adds a weighted-least-squares smoothing pass (cv2.ximgproc.createDisparityWLSFilter in opencv-contrib) that fills the invalidated holes from confident neighbors, guided by the image edges, a direct cousin of the edge-aware filtering in Chapter 7.

4. Reading a Disparity Map Like a Professional Intermediate

A disparity map is the most diagnostic image you will ever stare at: every artifact in it names its own cause, and learning to read those tells is the difference between trusting a depth sensor and being surprised by it on the road. Visualize disp with a colormap and inspect it against the failure catalog. Smooth surfaces should show smooth gradients; staircase banding indicates integer disparities (check the divide-by-16 and subpixel settings). A halo of invalid pixels hugging the left side of every foreground object is normal: those are right-camera occlusions, the geometry's fault, not the matcher's. Fat foreground edges that spill onto the background by half a block width are the edge-fattening signature; shrink blockSize or add the WLS filter. And a hard band of garbage along the left image border is the search range itself: a pixel at column $x < d_{\max}$ cannot have its match in view, so the leftmost numDisparities columns are structurally unmatchable. Production rigs simply crop them.