Depth from Two Views

This section was written by Ruimeng Yang (pull request #58), with help from an AI coding agent on the code. It reproduces the ideas of Chapter 40 of Foundations of Computer Vision by Antonio Torralba, Phillip Isola, and William T. Freeman. The book’s own figures are not reproduced here, because the book’s license covers only the work in full; links point to them instead. This section follows the order of the book’s stereo chapter and turns it into executable experiments. Where the book uses a photograph or a diagram, a link points to it. Where geometry, matching, failure analysis, or evaluation can be computed, the code computes it on a synthetic stereo pair with known ground truth, so every result can be measured. Each topic follows the same steps: the concept, the relevant figure in the book, the math, an experiment, its visualization, a quantitative reading of the result, and the failure cases.

import json
import math
import re
import time
from pathlib import Path

import matplotlib.colors as mcolors
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
import torch
import torch.nn.functional as F
from PIL import Image, ImageDraw

SEED = 7
_ = torch.manual_seed(SEED)
np.random.seed(SEED)
DEVICE = torch.device("cpu")

FOCAL_LENGTH_PX = 72.0
BASELINE_UNITS = 0.24
HEIGHT = 96
WIDTH = 144
DEFAULT_PATCH_SIZE = 7
DEFAULT_MAX_DISPARITY = 12
BAD_PIXEL_THRESHOLD = 1.0
EDGE_COLOR = "#0f766e"
ACCENT_RED = "#b91c1c"
ACCENT_BLUE = "#2563eb"
ACCENT_GOLD = "#d97706"

def benchmark_configuration(reference, target, patch_size, max_disparity, direction="left_to_right", warmups=2, measured=10):
    result = None
    for _ in range(warmups):
        result = run_match(reference, target, max_disparity=max_disparity, patch_size=patch_size, direction=direction)
    runtimes = []
    for _ in range(measured):
        start = time.perf_counter()
        result = run_match(reference, target, max_disparity=max_disparity, patch_size=patch_size, direction=direction)
        runtimes.append(time.perf_counter() - start)
    assert result is not None
    median_runtime = float(np.median(runtimes))
    return result, median_runtime, {
        "median_runtime_seconds": median_runtime,
        "min_runtime_seconds": float(np.min(runtimes)),
        "max_runtime_seconds": float(np.max(runtimes)),
        "warmup_runs": warmups,
        "measured_runs": measured,
    }

def as_float_tensor(array):
    if isinstance(array, torch.Tensor):
        return array.to(device=DEVICE, dtype=torch.float32)
    return torch.as_tensor(array, dtype=torch.float32, device=DEVICE)

def tensor_to_numpy(array):
    if isinstance(array, torch.Tensor):
        return array.detach().cpu().numpy()
    return np.asarray(array)

def masked_invalid_numpy(array):
    return np.ma.masked_invalid(tensor_to_numpy(array))

def finite_quantile(array, q):
    tensor = as_float_tensor(array)
    values = tensor[torch.isfinite(tensor)]
    return float(torch.quantile(values, q).item())

def finite_min(array):
    tensor = as_float_tensor(array)
    values = tensor[torch.isfinite(tensor)]
    return float(values.min().item())

def finite_max(array):
    tensor = as_float_tensor(array)
    values = tensor[torch.isfinite(tensor)]
    return float(values.max().item())

def box_filter(image, radius):
    image_t = as_float_tensor(image)
    if radius <= 0:
        return image_t.clone()
    kernel = 2 * radius + 1
    padded = F.pad(image_t[None, None], (radius, radius, radius, radius), mode="replicate")
    filtered = F.avg_pool2d(padded, kernel_size=kernel, stride=1)
    return filtered[0, 0]

def disparity_to_depth(disparity, focal_length=FOCAL_LENGTH_PX, baseline=BASELINE_UNITS):
    disparity_t = as_float_tensor(disparity)
    depth = torch.full_like(disparity_t, float("nan"))
    valid = torch.isfinite(disparity_t)
    depth[valid] = focal_length * baseline / torch.clamp(disparity_t[valid], min=1e-3)
    return depth if isinstance(disparity, torch.Tensor) else tensor_to_numpy(depth)

def depth_to_disparity(depth, focal_length=FOCAL_LENGTH_PX, baseline=BASELINE_UNITS):
    depth_t = as_float_tensor(depth)
    disparity = torch.full_like(depth_t, float("nan"))
    valid = depth_t > 1e-6
    disparity[valid] = focal_length * baseline / depth_t[valid]
    return disparity if isinstance(depth, torch.Tensor) else tensor_to_numpy(disparity)

def create_synthetic_scene(height=HEIGHT, width=WIDTH, focal_length=FOCAL_LENGTH_PX, baseline=BASELINE_UNITS):
    y, x = torch.meshgrid(
        torch.arange(height, dtype=torch.float32, device=DEVICE),
        torch.arange(width, dtype=torch.float32, device=DEVICE),
        indexing="ij",
    )
    left = 0.28 + 0.10 * torch.sin(x / 6.0) + 0.08 * torch.cos(y / 8.0) + 0.04 * torch.sin((x + 1.3 * y) / 11.0)
    disparity = torch.full((height, width), 2.2, dtype=torch.float32, device=DEVICE)
    region = np.full((height, width), "background", dtype=object)

    checker = (x > 16) & (x < 60) & (y > 18) & (y < 66)
    left[checker] = 0.38 + 0.24 * torch.remainder(torch.floor(x[checker] / 5.0) + torch.floor(y[checker] / 5.0), 2.0)
    disparity[checker] = 8.7
    region[tensor_to_numpy(checker)] = "checker"

    textureless = (x > 18) & (x < 64) & (y > 72) & (y < 91)
    left[textureless] = 0.63
    disparity[textureless] = 3.1
    region[tensor_to_numpy(textureless)] = "textureless"

    repeated = (x > 96) & (x < 138) & (y > 16) & (y < 52)
    left[repeated] = 0.34 + 0.22 * torch.remainder(torch.floor((x[repeated] - 96.0) / 4.0), 2.0)
    disparity[repeated] = 8.0
    region[tensor_to_numpy(repeated)] = "repeated"

    occluder = (x > 66) & (x < 76) & (y > 14) & (y < 82)
    left[occluder] = 0.86 - 0.10 * torch.remainder(torch.floor((y[occluder] - 14.0) / 6.0), 2.0)
    disparity[occluder] = 9.6
    region[tensor_to_numpy(occluder)] = "occluder"

    circle = (x - 88.0) ** 2 + (y - 38.0) ** 2 < 17.0 ** 2
    left[circle] = 0.22 + 0.55 * torch.exp(-((x[circle] - 88.0) ** 2 + (y[circle] - 38.0) ** 2) / 140.0)
    disparity[circle] = 4.6
    region[tensor_to_numpy(circle)] = "circle"

    ramp = (x > 78) & (x < 132) & (y > 60) & (y < 89)
    left[ramp] = 0.25 + 0.30 * ((x[ramp] - 78.0) / (132.0 - 78.0))
    disparity[ramp] = 3.5 + 2.8 * ((x[ramp] - 78.0) / (132.0 - 78.0))
    region[tensor_to_numpy(ramp)] = "ramp"

    left = torch.clamp(left, 0.0, 1.0)
    depth = disparity_to_depth(disparity, focal_length, baseline)

    right = torch.full_like(left, float("nan"))
    right_disp = torch.full_like(disparity, float("-inf"))
    source_visible = torch.zeros((height, width), dtype=torch.bool, device=DEVICE)
    invalid_projection = torch.zeros((height, width), dtype=torch.bool, device=DEVICE)
    source_xr = torch.full((height, width), -1, dtype=torch.long, device=DEVICE)

    for yy in range(height):
        for xx in range(width):
            d = float(disparity[yy, xx].item())
            xr = int(round(xx - d))
            source_xr[yy, xx] = xr
            if xr < 0 or xr >= width:
                invalid_projection[yy, xx] = True
                continue
            if d > float(right_disp[yy, xr].item()):
                right_disp[yy, xr] = d
                right[yy, xr] = left[yy, xx]

    row_x = torch.arange(width, dtype=torch.float32, device=DEVICE)
    for yy in range(height):
        valid_cols = torch.isfinite(right[yy])
        if int(valid_cols.sum().item()) >= 2:
            xp = row_x[valid_cols]
            fp = right[yy, valid_cols]
            right_idx = torch.searchsorted(xp, row_x)
            left_idx = torch.clamp(right_idx - 1, 0, xp.numel() - 1)
            right_idx = torch.clamp(right_idx, 0, xp.numel() - 1)
            x0 = xp[left_idx]
            x1 = xp[right_idx]
            y0 = fp[left_idx]
            y1 = fp[right_idx]
            denom = torch.where(torch.abs(x1 - x0) < 1e-6, torch.ones_like(x1), x1 - x0)
            alpha = torch.where(
                torch.abs(x1 - x0) < 1e-6,
                torch.zeros_like(row_x),
                (row_x - x0) / denom,
            )
            fill = y0 + alpha * (y1 - y0)
            right[yy] = fill
        elif int(valid_cols.sum().item()) == 1:
            right[yy] = right[yy, valid_cols][0]
        else:
            right[yy] = 0.0

    for yy in range(height):
        for xx in range(width):
            xr = int(source_xr[yy, xx].item())
            if xr < 0 or xr >= width:
                continue
            source_visible[yy, xx] = abs(float(right_disp[yy, xr].item()) - float(disparity[yy, xx].item())) < 1e-6

    occlusion_mask = ~source_visible
    return {
        "left": left,
        "right": right,
        "gt_disparity": disparity,
        "gt_depth": depth,
        "visible_mask": source_visible,
        "occlusion_mask": occlusion_mask,
        "invalid_projection": invalid_projection,
        "region": region,
    }

def shift_for_disparity(image, disparity, direction):
    image_t = as_float_tensor(image)
    shifted = torch.full_like(image_t, float("nan"))
    valid = torch.zeros_like(image_t, dtype=torch.bool)
    if disparity == 0:
        shifted[:] = image_t
        valid[:] = True
        return shifted, valid
    if direction == "left_to_right":
        shifted[:, disparity:] = image_t[:, :-disparity]
        valid[:, disparity:] = True
    elif direction == "right_to_left":
        shifted[:, :-disparity] = image_t[:, disparity:]
        valid[:, :-disparity] = True
    else:
        raise ValueError(f"Unsupported direction: {direction}")
    return shifted, valid

def build_support_mask(height, width, radius, max_disparity, direction):
    mask = torch.zeros((height, width), dtype=torch.bool, device=DEVICE)
    y0, y1 = radius, height - radius
    if direction == "left_to_right":
        x0, x1 = radius + max_disparity, width - radius
    elif direction == "right_to_left":
        x0, x1 = radius, width - radius - max_disparity
    else:
        raise ValueError(f"Unsupported direction: {direction}")
    if y1 > y0 and x1 > x0:
        mask[y0:y1, x0:x1] = True
    return mask

def build_cost_volume(reference, target, patch_size, max_disparity, direction):
    reference_t = as_float_tensor(reference)
    target_t = as_float_tensor(target)
    radius = patch_size // 2
    h, w = reference_t.shape
    cost_volume = torch.full((max_disparity + 1, h, w), float("nan"), dtype=torch.float32, device=DEVICE)
    support_mask = build_support_mask(h, w, radius, max_disparity, direction)
    for d in range(max_disparity + 1):
        aligned, aligned_valid = shift_for_disparity(target_t, d, direction)
        diff = (reference_t - aligned) ** 2
        diff = torch.where(aligned_valid, diff, torch.full_like(diff, 1e3))
        cost = box_filter(diff, radius)
        valid = support_mask & aligned_valid
        layer = torch.full((h, w), float("nan"), dtype=torch.float32, device=DEVICE)
        layer[valid] = cost[valid]
        cost_volume[d] = layer
    return cost_volume, support_mask

def disparity_from_cost_volume(cost_volume):
    cost_volume_t = as_float_tensor(cost_volume)
    filled = torch.where(torch.isfinite(cost_volume_t), cost_volume_t, torch.full_like(cost_volume_t, float("inf")))
    best = torch.argmin(filled, dim=0).to(torch.float32)
    valid = torch.isfinite(cost_volume_t).any(dim=0)
    best = torch.where(valid, best, torch.full_like(best, float("nan")))
    best_cost = torch.min(filled, dim=0).values
    best_cost = torch.where(valid, best_cost, torch.full_like(best_cost, float("nan")))
    return best, valid, best_cost

def compute_error_maps(pred_disparity, gt_disparity, gt_depth):
    pred_disparity_t = as_float_tensor(pred_disparity)
    gt_disparity_t = as_float_tensor(gt_disparity)
    gt_depth_t = as_float_tensor(gt_depth)
    disparity_error = torch.abs(pred_disparity_t - gt_disparity_t)
    pred_depth = as_float_tensor(disparity_to_depth(pred_disparity_t))
    depth_error = torch.abs(pred_depth - gt_depth_t)
    return disparity_error, pred_depth, depth_error

def summarize_metrics(pred_disparity, gt_disparity, gt_depth, eval_mask, runtime_seconds, consistency_mask=None):
    eval_mask_t = torch.as_tensor(eval_mask, dtype=torch.bool, device=DEVICE)
    disparity_error, pred_depth, depth_error = compute_error_maps(pred_disparity, gt_disparity, gt_depth)
    disp_values = disparity_error[eval_mask_t]
    useful_depth = eval_mask_t & torch.isfinite(depth_error) & torch.isfinite(as_float_tensor(pred_disparity)) & (as_float_tensor(pred_disparity) > 0.25)
    depth_values = depth_error[useful_depth]
    metrics = {
        "disparity_mae_px": float(disp_values.mean().item()),
        "bad_pixel_rate_gt_1px": float((disp_values > BAD_PIXEL_THRESHOLD).to(torch.float32).mean().item()),
        "valid_pixel_ratio": float(eval_mask_t.to(torch.float32).mean().item()),
        "depth_rmse_scene_units": float(torch.sqrt((depth_values ** 2).mean()).item()),
        "runtime_seconds": float(runtime_seconds),
    }
    if consistency_mask is not None:
        consistency_mask_t = torch.as_tensor(consistency_mask, dtype=torch.bool, device=DEVICE)
        metrics["left_right_consistency_rate"] = float(consistency_mask_t[eval_mask_t].to(torch.float32).mean().item())
    return metrics, disparity_error, pred_depth, depth_error

def run_match(reference, target, max_disparity, patch_size, direction="left_to_right"):
    start = time.perf_counter()
    cost_volume, support_mask = build_cost_volume(reference, target, patch_size, max_disparity, direction)
    pred_disparity, valid_mask, min_cost = disparity_from_cost_volume(cost_volume)
    runtime = time.perf_counter() - start
    return {
        "cost_volume": cost_volume,
        "support_mask": support_mask,
        "pred_disparity": pred_disparity,
        "valid_mask": valid_mask,
        "min_cost": min_cost,
        "runtime_seconds": runtime,
        "direction": direction,
    }

def sample_cost_curve(cost_volume, yy, xx):
    return as_float_tensor(cost_volume)[:, yy, xx]

def bilateral_brightness_variant(left, right):
    left_t = as_float_tensor(left)
    right_t = as_float_tensor(right)
    right_bright = torch.clamp(0.12 + 0.82 * right_t, 0.0, 1.0)
    return left_t.clone(), right_bright

def compute_right_to_left_consistency(left_disp, right_disp, left_valid_mask, right_valid_mask, eval_mask, tol=1.0):
    left_disp_t = as_float_tensor(left_disp)
    right_disp_t = as_float_tensor(right_disp)
    left_valid_t = torch.as_tensor(left_valid_mask, dtype=torch.bool, device=DEVICE)
    right_valid_t = torch.as_tensor(right_valid_mask, dtype=torch.bool, device=DEVICE)
    eval_mask_t = torch.as_tensor(eval_mask, dtype=torch.bool, device=DEVICE)
    h, w = left_disp_t.shape
    x_coords = torch.arange(w, device=DEVICE).view(1, w).expand(h, w)
    disp_indices = torch.round(torch.where(torch.isfinite(left_disp_t), left_disp_t, torch.zeros_like(left_disp_t))).to(torch.long)
    xr = x_coords - disp_indices
    in_bounds = (xr >= 0) & (xr < w)
    xr_safe = xr.clamp(0, w - 1)
    sampled_right_disp = torch.gather(right_disp_t, 1, xr_safe)
    sampled_right_valid = torch.gather(right_valid_t.to(torch.int64), 1, xr_safe).to(torch.bool)
    consistency = (
        eval_mask_t
        & left_valid_t
        & torch.isfinite(left_disp_t)
        & in_bounds
        & sampled_right_valid
        & torch.isfinite(sampled_right_disp)
        & (torch.abs(left_disp_t - sampled_right_disp) <= tol)
    )
    return consistency

def fit_subpixel_quadratic(cost_curve, best_disp):
    curve_t = as_float_tensor(cost_curve)
    if best_disp <= 0 or best_disp >= len(curve_t) - 1:
        return float(best_disp), None
    c1 = float(curve_t[best_disp - 1].item())
    c2 = float(curve_t[best_disp].item())
    c3 = float(curve_t[best_disp + 1].item())
    if not np.all(np.isfinite([c1, c2, c3])):
        return float(best_disp), None
    denom = c1 - 2.0 * c2 + c3
    if abs(denom) < 1e-9:
        return float(best_disp), None
    offset = 0.5 * (c1 - c3) / denom
    refined = float(best_disp + offset)
    xs = np.linspace(best_disp - 1, best_disp + 1, 200)
    ys = c2 + 0.5 * denom * (xs - best_disp) ** 2 + 0.5 * (c3 - c1) * (xs - best_disp)
    return refined, (xs, ys)

def create_constant_disparity_pair(height=72, width=112, disparity_px=4):
    y, x = torch.meshgrid(
        torch.arange(height, dtype=torch.float32, device=DEVICE),
        torch.arange(width, dtype=torch.float32, device=DEVICE),
        indexing="ij",
    )
    left = 0.28 + 0.22 * torch.sin(x / 3.7) + 0.18 * torch.cos(y / 5.1) + 0.11 * torch.sin((1.9 * x + 0.8 * y) / 6.3)
    left = torch.clamp(left, 0.0, 1.0)
    right = torch.zeros_like(left)
    right[:, :-disparity_px] = left[:, disparity_px:]
    left_visible_mask = torch.zeros((height, width), dtype=torch.bool, device=DEVICE)
    left_visible_mask[:, disparity_px:] = True
    right_visible_mask = torch.zeros((height, width), dtype=torch.bool, device=DEVICE)
    right_visible_mask[:, : width - disparity_px] = True
    return {
        "left": left,
        "right": right,
        "left_visible_mask": left_visible_mask,
        "right_visible_mask": right_visible_mask,
        "disparity_px": float(disparity_px),
    }

def run_constant_disparity_sanity_test():
    sanity = create_constant_disparity_pair()
    left_match = run_match(
        sanity["left"],
        sanity["right"],
        max_disparity=8,
        patch_size=7,
        direction="left_to_right",
    )
    right_match = run_match(
        sanity["right"],
        sanity["left"],
        max_disparity=8,
        patch_size=7,
        direction="right_to_left",
    )
    left_eval_mask = sanity["left_visible_mask"] & left_match["valid_mask"]
    right_eval_mask = sanity["right_visible_mask"] & right_match["valid_mask"]
    consistency = compute_right_to_left_consistency(
        left_match["pred_disparity"],
        right_match["pred_disparity"],
        left_match["valid_mask"],
        right_match["valid_mask"],
        left_eval_mask,
        tol=0.5,
    )
    median_left = float(torch.median(left_match["pred_disparity"][left_eval_mask]).item())
    median_right = float(torch.median(right_match["pred_disparity"][right_eval_mask]).item())
    consistency_rate = float(consistency[left_eval_mask].to(torch.float32).mean().item())
    print(
        f"Controlled 4 px sanity test: median left-to-right disparity = {median_left:.3f} px, "
        f"median right-to-left disparity = {median_right:.3f} px, "
        f"left-right consistency rate = {consistency_rate:.3f}"
    )
    assert abs(median_left - 4.0) <= 0.25
    assert abs(median_right - 4.0) <= 0.25
    assert consistency_rate >= 0.95
    return {
        "median_left_disparity": median_left,
        "median_right_disparity": median_right,
        "left_right_consistency_rate": consistency_rate,
    }

def warp_left_to_right_fractional(left_image, disparity_px):
    left_t = as_float_tensor(left_image)
    h, w = left_t.shape
    ys, xs = torch.meshgrid(
        torch.arange(h, dtype=torch.float32, device=DEVICE),
        torch.arange(w, dtype=torch.float32, device=DEVICE),
        indexing="ij",
    )
    source_x = xs + disparity_px
    source_y = ys
    valid = (source_x >= 0.0) & (source_x <= (w - 1))
    grid_x = 2.0 * source_x / max(w - 1, 1) - 1.0
    grid_y = 2.0 * source_y / max(h - 1, 1) - 1.0
    grid = torch.stack((grid_x, grid_y), dim=-1)[None]
    warped = F.grid_sample(
        left_t[None, None],
        grid,
        mode="bilinear",
        padding_mode="zeros",
        align_corners=True,
    )[0, 0]
    return warped, valid

def run_fractional_disparity_experiment(ground_truth_disparity=8.7):
    height, width = 48, 96
    y, x = torch.meshgrid(
        torch.arange(height, dtype=torch.float32, device=DEVICE),
        torch.arange(width, dtype=torch.float32, device=DEVICE),
        indexing="ij",
    )
    left = 0.31 + 0.24 * torch.sin(x / 2.9) + 0.19 * torch.cos(y / 4.1) + 0.12 * torch.sin((1.4 * x + 0.7 * y) / 5.3)
    left = torch.clamp(left, 0.0, 1.0)
    right, valid = warp_left_to_right_fractional(left, ground_truth_disparity)
    match = run_match(left, right, max_disparity=12, patch_size=9, direction="left_to_right")
    yy, xx = 22, 48
    assert bool(valid[yy, xx].item()), "Fractional-disparity sample point fell outside the valid warp support."
    curve = sample_cost_curve(match["cost_volume"], yy, xx)
    best_disp = int(torch.argmin(torch.where(torch.isfinite(curve), curve, torch.full_like(curve, float("inf")))).item())
    refined_disp, fitted = fit_subpixel_quadratic(curve, best_disp)
    integer_error = abs(best_disp - ground_truth_disparity)
    refined_error = abs(refined_disp - ground_truth_disparity)
    print(
        f"Fractional disparity experiment: gt = {ground_truth_disparity:.3f} px, "
        f"integer estimate = {best_disp:.3f} px, refined estimate = {refined_disp:.3f} px, "
        f"integer error = {integer_error:.3f} px, refined error = {refined_error:.3f} px"
    )
    assert refined_error < integer_error
    assert refined_error < 0.15
    return {
        "left": left,
        "right": right,
        "valid_mask": valid,
        "curve": curve,
        "yy": yy,
        "xx": xx,
        "ground_truth_disparity": float(ground_truth_disparity),
        "integer_disparity": float(best_disp),
        "refined_disparity": float(refined_disp),
        "integer_error": float(integer_error),
        "refined_error": float(refined_error),
        "fitted_curve": fitted,
    }

def rot_x(theta):
    c, s = np.cos(theta), np.sin(theta)
    return np.array([[1.0, 0.0, 0.0], [0.0, c, -s], [0.0, s, c]], dtype=np.float64)

def rot_y(theta):
    c, s = np.cos(theta), np.sin(theta)
    return np.array([[c, 0.0, s], [0.0, 1.0, 0.0], [-s, 0.0, c]], dtype=np.float64)

def skew(vec):
    tx, ty, tz = vec
    return np.array([[0.0, -tz, ty], [tz, 0.0, -tx], [-ty, tx, 0.0]], dtype=np.float64)

def project_points(K, R, t, points):
    pixels = []
    for point in points:
        cam = R @ point + t
        pix = K @ cam
        pixels.append(np.array([pix[0] / pix[2], pix[1] / pix[2], 1.0], dtype=np.float64))
    return np.stack(pixels)

def line_segment_in_frame(line, width, height):
    a, b, c = line
    points = []
    for x in [0.0, width - 1.0]:
        if abs(b) > 1e-9:
            y = -(a * x + c) / b
            if -10.0 <= y <= height + 10.0:
                points.append((x, y))
    for y in [0.0, height - 1.0]:
        if abs(a) > 1e-9:
            x = -(b * y + c) / a
            if -10.0 <= x <= width + 10.0:
                points.append((x, y))
    unique = []
    for point in points:
        if all((point[0] - other[0]) ** 2 + (point[1] - other[1]) ** 2 > 1e-6 for other in unique):
            unique.append(point)
    return unique[:2] if len(unique) >= 2 else None

scene = create_synthetic_scene()
eval_mask_template = scene["visible_mask"]
global_artifacts = {"sanity_test": run_constant_disparity_sanity_test()}

Controlled 4 px sanity test: median left-to-right disparity = 4.000 px, median right-to-left disparity = 4.000 px, left-right consistency rate = 0.959

from IPython.display import display


def show_figure(fig):
    """Render a finished figure inline and release it."""
    display(fig)
    plt.close(fig)

Introduction

Stereo begins as a perceptual phenomenon: the left and right eyes see slightly different image positions, and those displacements can be interpreted as depth. Figure 40.1 in the book: stereo anaglyph of the Titanic and the red/cyan viewing setup that turns left-right displacement into a 3D percept. The chapter quickly turns that perception into a two-part computational problem:

geometry: where can a corresponding point lie?
matching: which candidate point is the correct one?

The rest of this section follows the book’s order and adds executable experiments for the same ideas.

Stereo cues

How far away is a boat?

The boat example introduces two depth cues before stereo algorithms appear:

a single-view estimate using observer height h and angle alpha, with d = h / tan(alpha);
a two-view triangulation estimate using baseline t and angles alpha and beta, with d = t sin(alpha) sin(beta) / sin(alpha + beta).

Figure 40.2 in the book: two boat-distance constructions: one using a horizon reference and one using triangulation between two observation points.

Depth from image disparities

The random-dot stereogram isolates disparity from recognizable object identity. It shows that a depth percept can arise purely from left-right displacement. Figure 40.3 in the book: random-dot stereogram showing that disparity alone can create a depth percept.

Building a stereo pinhole camera

The chapter also grounds the perceptual story in image formation by showing a homemade anaglyph pinhole camera. Figure 40.4 in the book: an anaglyph pinhole camera built from two pinholes, color filters, and a projection plane. The next cell turns those ideas into five figures: geometric stereo cues, rectified stereo geometry, the inverse disparity-depth curve, far-depth sensitivity, and baseline sensitivity.

single_height = 30.0
single_alpha_deg = 8.0
single_distance = single_height / math.tan(math.radians(single_alpha_deg))

baseline_t = 120.0
alpha_deg = 38.0
beta_deg = 29.0
triangulated_distance = baseline_t * math.sin(math.radians(alpha_deg)) * math.sin(math.radians(beta_deg)) / math.sin(
    math.radians(alpha_deg + beta_deg)
)

def make_random_dot_stereogram(height=84, width=128, square_size=28, shift=7):
    base = (np.random.rand(height, width) > 0.5).astype(np.float32)
    left = base.copy()
    right = base.copy()
    y0 = height // 2 - square_size // 2
    x0 = width // 2 - square_size // 2
    patch = (np.random.rand(square_size, square_size) > 0.5).astype(np.float32)
    left[y0 : y0 + square_size, x0 : x0 + square_size] = patch
    right[y0 : y0 + square_size, x0 - shift : x0 - shift + square_size] = patch
    return left, right, shift

stereogram_left, stereogram_right, stereogram_shift = make_random_dot_stereogram()

disparities = np.linspace(0.5, 18.0, 300)
depths = FOCAL_LENGTH_PX * BASELINE_UNITS / disparities

depth_grid = np.linspace(1.5, 11.0, 240)
disparity_nominal = depth_to_disparity(depth_grid)
disparity_errors = [0.25, 0.5, 1.0]

baselines = [0.12, 0.24, 0.36]

global_artifacts["geometry_summary"] = {
    "single_view_distance": single_distance,
    "triangulated_distance": triangulated_distance,
    "stereogram_shift_px": int(stereogram_shift),
}
print(json.dumps(global_artifacts["geometry_summary"], indent=2))

{
  "single_view_distance": 213.46109167152625,
  "triangulated_distance": 38.91063973428839,
  "stereogram_shift_px": 7
}

The figures above show, in order:

Boat triangulation and a random-dot stereogram make the book’s depth cue visible before any stereo algorithm appears.
Rectified stereo geometry with baseline, focal length, corresponding points, disparity, and a triangulated 3D point.
The inverse relationship between disparity and depth. Equal disparity steps do not correspond to equal depth steps.
A fixed disparity error creates much larger depth error at long range, which is why far geometry is fragile.
Larger baselines improve disparity signal but also create stronger view changes, overlap loss, and occlusion risk.

These figures sharpen three points from the book:

positive disparity means the corresponding point appears farther left in the right image;
depth is nonlinear in disparity;
far points are especially sensitive to small disparity errors, so geometry alone does not make stereo easy.

Model-based methods

Triangulation

The simple rectified geometry above is the special case used to derive the book’s core equations:

disparity: d = x_L - x_R
depth: Z = fB / d

Figure 40.5 in the book: rectified stereo geometry with focal length, baseline, image coordinates, and a triangulated 3D point.

Stereo matching

Once the geometry is fixed, the practical problem is correspondence. Which pixel in the right image matches a given pixel in the left? Figure 40.6 in the book: a real stereo office pair with corresponding features and their displacements highlighted. The next experiments use a synthetic rectified pair to compute a cost volume, run winner-takes-all disparity estimation, and measure actual error.

scene = create_synthetic_scene()
left = scene["left"]
right = scene["right"]
gt_disparity = scene["gt_disparity"]
gt_depth = scene["gt_depth"]
visible_mask = scene["visible_mask"]
region = scene["region"]

sample_y, sample_x = 36, 44
pixel_match = run_match(left, right, patch_size=1, max_disparity=DEFAULT_MAX_DISPARITY, direction="left_to_right")
patch_match = run_match(left, right, patch_size=DEFAULT_PATCH_SIZE, max_disparity=DEFAULT_MAX_DISPARITY, direction="left_to_right")
bright_left, bright_right = bilateral_brightness_variant(left, right)
bright_match = run_match(bright_left, bright_right, patch_size=DEFAULT_PATCH_SIZE, max_disparity=DEFAULT_MAX_DISPARITY, direction="left_to_right")

disparity_error, pred_depth, depth_error = compute_error_maps(
    patch_match["pred_disparity"], gt_disparity, gt_depth
)
eval_mask = visible_mask & patch_match["valid_mask"]
metrics_base, disparity_error, pred_depth, depth_error = summarize_metrics(
    patch_match["pred_disparity"], gt_disparity, gt_depth, eval_mask, patch_match["runtime_seconds"]
)

patch_sizes = [3, 5, 7, 9]
patch_results = []
for patch_size in patch_sizes:
    result = run_match(left, right, patch_size=patch_size, max_disparity=DEFAULT_MAX_DISPARITY, direction="left_to_right")
    mask = visible_mask & result["valid_mask"]
    metrics, _, _, _ = summarize_metrics(result["pred_disparity"], gt_disparity, gt_depth, mask, result["runtime_seconds"])
    patch_results.append((patch_size, result, metrics))

max_ranges = [6, 10, 14, 18]
range_results = []

global_artifacts["base_match"] = patch_match
global_artifacts["patch_results"] = patch_results
global_artifacts["range_results"] = range_results
global_artifacts["metrics_base"] = metrics_base
global_artifacts["scene_summary"] = {
    "ground_truth_disparity_range_px": [finite_min(gt_disparity), finite_max(gt_disparity)],
    "visible_pixel_count": int(visible_mask.to(torch.int64).sum().item()),
    "base_metrics": metrics_base,
}
print(json.dumps(global_artifacts["scene_summary"], indent=2))

{
  "ground_truth_disparity_range_px": [
    2.200000047683716,
    9.600000381469727
  ],
  "visible_pixel_count": 12596,
  "base_metrics": {
    "disparity_mae_px": 1.4484385251998901,
    "bad_pixel_rate_gt_1px": 0.3103710114955902,
    "valid_pixel_ratio": 0.7623698115348816,
    "depth_rmse_scene_units": 1.6601009368896484,
    "runtime_seconds": 0.005328660976374522
  }
}

The figures above show, in order:

Single-pixel matching is brittle; patch aggregation stabilizes the cost surface, while brightness shifts still move the optimum.
Each disparity slice of the cost volume answers one question: how plausible is this disparity at every image location?
Winner-takes-all disparity estimation chooses d*(x,y) = argmin_d C(x,y,d) independently at each pixel.
Patch size is a real tradeoff: too small is noisy, too large blurs across discontinuities and occlusions.
The disparity search range must be large enough to include valid solutions but small enough to avoid wasted runtime and extra ambiguity.

Quantitatively, the experiments expose the actual matching objective: d*(x, y) = argmin_d C(x, y, d) where C(x, y, d) is the matching cost at pixel (x, y) for candidate disparity d. The sweeps make three failure modes visible:

too-small patches are unstable;
too-large patches bleed across depth discontinuities;
incorrect disparity ranges either truncate the solution or waste compute.

Finding image features

Feature-based stereo is the book’s answer to the fragility of raw intensity matching. Figure 40.8 in the book: feature detections on the office pair and the interpolated depth result. A good feature is localizable under small translations. The chapter expresses that with the Harris patch energy E(Delta x, Delta y) = sum_(x,y in P) (l(x,y) - l(x + Delta x, y + Delta y))^2 and then motivates richer local descriptors.

Local image descriptors

Figure 40.9 in the book: orientation-based local descriptors and the spatial pooling idea behind SIFT-style matching. Oriented local structure is often more stable than raw intensity values, especially under small view changes.

Interpolation between feature matches

Even after sparse feature matching, a system still needs interpolation or regularization to obtain dense depth. The next figures analyze the failure cases that make that interpolation necessary.

base_match = global_artifacts["base_match"]
pred_left = base_match["pred_disparity"]
texture_y, texture_x = 80, 30
repeat_y, repeat_x = 30, 112
subpixel_y, subpixel_x = 36, 44

right_to_left = run_match(right, left, patch_size=DEFAULT_PATCH_SIZE, max_disparity=DEFAULT_MAX_DISPARITY, direction="right_to_left")
consistency_eval_mask = visible_mask & base_match["valid_mask"]
consistency = compute_right_to_left_consistency(
    pred_left,
    right_to_left["pred_disparity"],
    base_match["valid_mask"],
    right_to_left["valid_mask"],
    consistency_eval_mask,
    tol=1.0,
)
inconsistency = consistency_eval_mask & ~consistency

fractional_demo = run_fractional_disparity_experiment(ground_truth_disparity=8.7)
curve = fractional_demo["curve"]
best_disp = fractional_demo["integer_disparity"]
refined_disp = fractional_demo["refined_disparity"]
fitted = fractional_demo["fitted_curve"]

Fractional disparity experiment: gt = 8.700 px, integer estimate = 9.000 px, refined estimate = 8.691 px, integer error = 0.300 px, refined error = 0.009 px

global_artifacts["consistency_mask"] = consistency
global_artifacts["right_to_left"] = right_to_left
global_artifacts["subpixel_demo"] = {
    "integer_disparity": float(best_disp),
    "refined_disparity": float(refined_disp),
    "ground_truth_disparity": fractional_demo["ground_truth_disparity"],
    "integer_error": fractional_demo["integer_error"],
    "refined_error": fractional_demo["refined_error"],
}
print(json.dumps({"subpixel_demo": global_artifacts["subpixel_demo"]}, indent=2))

{
  "subpixel_demo": {
    "integer_disparity": 9.0,
    "refined_disparity": 8.690843069220339,
    "ground_truth_disparity": 8.7,
    "integer_error": 0.3000000000000007,
    "refined_error": 0.0091569307796604
  }
}

The figures above show, in order:

In a low-texture region the matching-cost surface becomes flat, so many disparities look almost equally plausible.
Repeated patterns create several near-identical alignments, which makes the cost surface multimodal.
Left-right consistency reveals pixels whose correspondences are unstable or absent because of occlusion.
A local parabola fit around the best integer disparity can recover a subpixel estimate when the cost curve is well behaved.

The failure modes, side by side:

textureless regions fail because the cost surface is flat;
repetitive textures fail because several disparities have similar cost;
occlusions fail because a point may not exist in both views at all;
brightness changes shift the entire cost curve, even when geometry is correct;
too-small patches are noisy;
too-large patches cross depth boundaries and mix objects;
far-depth sensitivity amplifies small disparity mistakes into large depth errors.

Practical mitigations exist, but each one leaves a tradeoff behind. Smoother costs reduce noise but blur discontinuities, larger baselines improve precision but increase occlusion, and subpixel fits help only when the local minimum is already reliable.

Constraints for arbitrary cameras

Once you leave the rectified special case, corresponding points are no longer found by searching the same row in the second image. Figure 40.10 in the book: candidate-match ambiguity for a point viewed under arbitrary camera geometry. Figure 40.11 in the book: the viewing ray from camera 1 projects to an epipolar line in camera 2. Figure 40.12 in the book: epipolar plane, epipolar lines, and epipoles for a stereo pair. The epipolar constraint is the algebraic form of that geometry: x'^T F x = 0 Valid correspondences should make the residual close to zero.

The essential and fundamental matrices

The essential matrix uses calibrated camera coordinates, while the fundamental matrix absorbs the intrinsic calibration and operates in image coordinates.

Epipolar lines: the game

Figure 40.13 in the book: epipolar-line intuition game matching camera arrangements to line families.

Image rectification

Rectification is a practical warp that makes epipolar lines horizontal again, reducing a 2D search problem to a 1D search problem.

width_ep, height_ep = 320, 240
K = np.array([[220.0, 0.0, 160.0], [0.0, 220.0, 120.0], [0.0, 0.0, 1.0]], dtype=np.float64)
R1 = np.eye(3)
t1 = np.zeros(3)
R2 = rot_y(np.deg2rad(12.0)) @ rot_x(np.deg2rad(-5.0))
t2 = np.array([0.55, 0.06, 0.02], dtype=np.float64)

world_points = np.array(
    [[0.00, 0.00, 3.2], [0.20, -0.05, 3.8], [-0.28, 0.08, 4.1], [0.18, 0.14, 2.9]],
    dtype=np.float64,
)
p1 = project_points(K, R1, t1, world_points)
p2 = project_points(K, R2, t2, world_points)
E = skew(t2) @ R2
Fmat = np.linalg.inv(K).T @ E @ np.linalg.inv(K)
residuals = [float(p2[idx] @ Fmat @ p1[idx]) for idx in range(len(world_points))]

global_artifacts["epipolar"] = {"residuals": residuals}
print(json.dumps(global_artifacts["epipolar"], indent=2))

{
  "residuals": [
    0.0,
    2.7755575615628914e-17,
    -5.551115123125783e-17,
    -5.551115123125783e-17
  ]
}

The figures above show, in order:

The point-to-line relation becomes algebraic through the epipolar constraint x'^T F x = 0.
Conceptual rectification illustration: it shows the geometric goal of horizontal scanlines, not the output of a full calibrated rectification pipeline.

These figures reconnect the later epipolar geometry back to the earlier rectified stereo experiments. Rectification is not a different problem; it is a practical reparameterization of the same correspondence constraint.

Learning-based methods

The book does not present a full deep-learning survey. Instead, it emphasizes why many learned systems still predict disparity rather than depth directly, and why they often use a rectified cost-volume pipeline. Figure 40.14 in the book: two-stage CNN stereo pipeline: feature extraction, cost volume, cost aggregation, and disparity estimate. The ideas to keep from it:

disparity is easier to regularize than depth under a fixed stereo rig;
cost volumes remain central even in learned stereo;
the network predicts disparity, but geometry still interprets the result.

Evaluation

The evaluation is executable: it measures classical block matching on the synthetic stereo pair.

base_match = global_artifacts["base_match"]
pred_disparity = base_match["pred_disparity"]
valid_mask = base_match["valid_mask"]
consistency_mask = global_artifacts["consistency_mask"]
eval_mask = visible_mask & valid_mask
metrics_final, disparity_error, pred_depth, depth_error = summarize_metrics(
    pred_disparity,
    gt_disparity,
    gt_depth,
    eval_mask,
    base_match["runtime_seconds"],
    consistency_mask=consistency_mask,
)

operating_points = []
for patch_size in [3, 5, 7, 9]:
    result, median_runtime, _ = benchmark_configuration(
        left,
        right,
        patch_size=patch_size,
        max_disparity=DEFAULT_MAX_DISPARITY,
        direction="left_to_right",
    )
    mask = visible_mask & result["valid_mask"]
    metrics, _, _, _ = summarize_metrics(result["pred_disparity"], gt_disparity, gt_depth, mask, median_runtime)
    operating_points.append(
        {
            "label": f"patch={patch_size}, range={DEFAULT_MAX_DISPARITY}",
            "runtime_seconds": float(metrics["runtime_seconds"]),
            "disparity_mae_px": float(metrics["disparity_mae_px"]),
            "family": "patch sweep",
        }
    )
for max_disp in [6, 10, 14, 18]:
    result, median_runtime, _ = benchmark_configuration(
        left,
        right,
        patch_size=DEFAULT_PATCH_SIZE,
        max_disparity=max_disp,
        direction="left_to_right",
    )
    mask = visible_mask & result["valid_mask"]
    metrics, _, _, _ = summarize_metrics(result["pred_disparity"], gt_disparity, gt_depth, mask, median_runtime)
    operating_points.append(
        {
            "label": f"patch={DEFAULT_PATCH_SIZE}, range={max_disp}",
            "runtime_seconds": float(metrics["runtime_seconds"]),
            "disparity_mae_px": float(metrics["disparity_mae_px"]),
            "family": "range sweep",
        }
    )
operating_points.sort(key=lambda row: (row["runtime_seconds"], row["disparity_mae_px"]))

global_artifacts["metrics_final"] = metrics_final
global_artifacts["operating_points"] = operating_points
print("Results table")
print("-" * 74)
print(f"{'label':<28} {'runtime(s)':>12} {'MAE(px)':>10}")
print("-" * 74)
for row in operating_points:
    print(f"{row['label']:<28} {row['runtime_seconds']:>12.3f} {row['disparity_mae_px']:>10.3f}")
print("-" * 74)
print(json.dumps({"metrics_final": metrics_final}, indent=2))

Results table
--------------------------------------------------------------------------
label                          runtime(s)    MAE(px)
--------------------------------------------------------------------------
patch=7, range=6                    0.003      2.992
patch=7, range=10                   0.004      1.424
patch=7, range=12                   0.005      1.448
patch=9, range=12                   0.005      1.585
patch=3, range=12                   0.005      1.775
patch=5, range=12                   0.005      1.346
patch=7, range=14                   0.005      1.482
patch=7, range=18                   0.006      1.589
--------------------------------------------------------------------------
{
  "metrics_final": {
    "disparity_mae_px": 1.4484385251998901,
    "bad_pixel_rate_gt_1px": 0.3103710114955902,
    "valid_pixel_ratio": 0.7623698115348816,
    "depth_rmse_scene_units": 1.6601009368896484,
    "runtime_seconds": 0.005328660976374522,
    "left_right_consistency_rate": 0.837365984916687
  }
}

The figures above show, in order:

Evaluation should expose both disparity-space and depth-space error, because small disparity errors can turn into large depth errors far from the cameras.
Runtime and accuracy should be read together. The timings are machine-dependent and are intended only for relative comparison within this run.

The code reports these metrics:

disparity MAE: mean absolute disparity error over valid visible pixels;
bad-pixel rate: fraction of valid pixels with absolute disparity error greater than 1 pixel;
valid-pixel ratio: fraction of all pixels that remain usable after support and visibility checks;
depth RMSE: root-mean-square error after converting disparity to depth;
left-right consistency rate: fraction of valid pixels that agree with a reverse-direction disparity check;
runtime: wall-clock runtime for the selected matching configuration.

Concluding remarks

The book ends by staying realistic about stereo:

correspondence is hard even in the rectified case;
far geometry is fragile because depth is nonlinear in disparity;
occlusions and repeated structure create genuine ambiguity;
practical systems mix geometry, matching, regularization, and evaluation rather than relying on a single elegant formula.

The experiments in this section reflect that balance: they run the book’s ideas as code and measure the consequences directly.

Edit this page on GitHub or file an issue.

Physical AI

3D Reconstruction

Imitation Learning

World Models

VLA Models

Sim-to-Real Transfer

Depth from Two Views

Introduction

Stereo cues

How far away is a boat?

Depth from image disparities

Building a stereo pinhole camera

Model-based methods

Triangulation

Stereo matching

Finding image features

Local image descriptors

Interpolation between feature matches

Constraints for arbitrary cameras

The essential and fundamental matrices

Epipolar lines: the game

Image rectification

Learning-based methods

Evaluation

Concluding remarks

​Introduction

​Stereo cues

​How far away is a boat?

​Depth from image disparities

​Building a stereo pinhole camera

​Model-based methods

​Triangulation

​Stereo matching

​Finding image features

​Local image descriptors

​Interpolation between feature matches

​Constraints for arbitrary cameras

​The essential and fundamental matrices

​Epipolar lines: the game

​Image rectification

​Learning-based methods

​Evaluation

​Concluding remarks

Introduction

Stereo cues

How far away is a boat?

Depth from image disparities

Building a stereo pinhole camera

Model-based methods

Triangulation

Stereo matching

Finding image features

Local image descriptors

Interpolation between feature matches

Constraints for arbitrary cameras

The essential and fundamental matrices

Epipolar lines: the game

Image rectification

Learning-based methods

Evaluation

Concluding remarks