{ "cells": [ { "cell_type": "markdown", "id": "d0f27f66", "metadata": {}, "source": [ "# Lesson 21: Stereo Matching\n", "\n", "Two cameras viewing the same scene from slightly different positions see the same 3D points shifted by different amounts depending on depth — nearby points shift more, distant points shift less. **Stereo matching** finds these shifts (the **disparity**) at every pixel, which convert into depth via triangulation. The key problem is correspondence: finding, for every pixel in the left image, the matching pixel in the right. This lesson builds a **block matcher** that solves it by comparing small patches with sum-of-squared differences (SSD) — a close cousin of Lesson 20's optical flow, but restricted to a 1D search along a single row instead of a full 2D neighborhood." ] }, { "cell_type": "code", "execution_count": null, "id": "98c8cf08", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import cv2\n", "import matplotlib.pyplot as plt\n", "\n", "# shared colormap: invalid/masked-out disparities (NaN) render as gray instead of viridis's default\n", "invalid_cmap = plt.cm.viridis.copy()\n", "invalid_cmap.set_bad('gray')" ] }, { "cell_type": "markdown", "id": "70662447", "metadata": {}, "source": [ "## Rectified stereo: why the search is 1D\n", "\n", "For a pair of cameras that are side-by-side, pointed the same direction, with parallel image planes (a **rectified** stereo pair — real camera rigs are calibrated and warped to approximate this), a fundamental fact of epipolar geometry applies: the corresponding point for any pixel in the left image lies **on the same row** in the right image. This collapses the search for a match from a 2D image search down to a 1D scan along one row, and the horizontal offset between the two matching positions is the **disparity** $d$.\n", "\n", "Disparity relates to depth by\n", "\n", "$$Z = \\frac{f \\cdot B}{d}$$\n", "\n", "where $f$ is the focal length and $B$ is the baseline (distance between the two camera centers). For a rectified pair, depth is **inversely proportional to disparity**: nearby objects have large disparity, distant objects have small disparity, and an object infinitely far away has zero disparity." ] }, { "cell_type": "markdown", "id": "582b1243", "metadata": {}, "source": [ "## A synthetic stereo pair with known ground truth\n", "\n", "We build a left image of pure random texture (so every patch is locally distinctive — no aperture-problem ambiguity) and construct the right image by shifting each pixel left by its true disparity, which we set to three different constant values for three depth \"planes\": a background and two nearer rectangles. Nearer surfaces shift more than farther ones, so they uncover background pixels on one side that have no corresponding pixel in the left image at all — an **occlusion**, which we patch with fresh random noise so every right-image pixel still has some value." ] }, { "cell_type": "code", "execution_count": null, "id": "d12a2602", "metadata": {}, "outputs": [], "source": [ "rng = np.random.default_rng(0)\n", "h, w = 150, 200\n", "left = rng.integers(0, 255, (h, w)).astype(np.uint8)\n", "\n", "true_disparity = np.full((h, w), 5, dtype=np.int32) # background (disparity = 5)\n", "true_disparity[30:120, 50:150] = 15 # a close rectangle (disparity = 15)\n", "true_disparity[60:90, 80:120] = 25 # a closer rectangle (disparity = 25)\n", "\n", "right = np.full((h, w), -1, dtype=np.int32)\n", "for y in range(h):\n", " for x in range(w):\n", " xr = x - true_disparity[y, x]\n", " if 0 <= xr < w:\n", " right[y, xr] = left[y, x]\n", "occluded = right == -1 # positions no left pixel maps to: disocclusions, filled with fresh noise\n", "right[occluded] = rng.integers(0, 255, occluded.sum())\n", "right = right.astype(np.uint8)\n", "\n", "fig, axes = plt.subplots(1, 3, figsize=(10, 3.5))\n", "axes[0].imshow(left, cmap='gray')\n", "axes[0].set_title('Left image', fontsize=9)\n", "axes[1].imshow(right, cmap='gray')\n", "axes[1].set_title('Right image', fontsize=9)\n", "im = axes[2].imshow(true_disparity, cmap='viridis')\n", "axes[2].set_title('Ground-truth disparity', fontsize=9)\n", "for ax in axes:\n", " ax.axis('off')\n", "plt.colorbar(im, ax=axes, orientation='vertical', fraction=0.046, pad=0.04, label='disparity (px)')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "a8ec13c0", "metadata": {}, "source": [ "## Block matching along the scanline\n", "\n", "For each pixel in the left image, we slide a small window along the *same row* of the right image over a range of candidate disparities and keep the disparity with the lowest sum-of-squared-differences. The implementation below has three nested Python loops, so it's slow — we run it on a small crop rather than the full image." ] }, { "cell_type": "code", "execution_count": null, "id": "527a664b", "metadata": {}, "outputs": [], "source": [ "def block_match_stereo(left, right, block_size=7, max_disp=30):\n", " h, w = left.shape\n", " half = block_size // 2\n", " left_f, right_f = left.astype(np.float64), right.astype(np.float64)\n", " disparity_map = np.zeros((h, w), dtype=np.float64)\n", "\n", " for y in range(half, h - half):\n", " for x in range(half, w - half):\n", " left_patch = left_f[y - half:y + half + 1, x - half:x + half + 1]\n", " best_d, best_cost = 0, np.inf\n", " for d in range(max_disp + 1):\n", " xr = x - d\n", " if xr - half < 0:\n", " break\n", " right_patch = right_f[y - half:y + half + 1, xr - half:xr + half + 1]\n", " cost = ((left_patch - right_patch)**2).sum()\n", " if cost < best_cost:\n", " best_cost, best_d = cost, d\n", " disparity_map[y, x] = best_d\n", " return disparity_map\n", "\n", "# a small crop, since the pure-Python pixel loop is slow\n", "crop = np.s_[0:100, 0:100]\n", "manual_disp = block_match_stereo(left[crop], right[crop], block_size=7, max_disp=30)\n", "\n", "fig, axes = plt.subplots(1, 3, figsize=(10, 3.5))\n", "axes[0].imshow(left[crop], cmap='gray')\n", "axes[0].set_title(f'Left crop ({left[crop].shape[0]}x{left[crop].shape[1]})', fontsize=9)\n", "axes[1].imshow(true_disparity[crop], cmap='viridis', vmin=0, vmax=30)\n", "axes[1].set_title(f'True disparity ({true_disparity[crop].shape[0]}x{true_disparity[crop].shape[1]})', fontsize=9)\n", "im = axes[2].imshow(manual_disp, cmap='viridis', vmin=0, vmax=30)\n", "axes[2].set_title(f'Recovered disparity ({manual_disp.shape[0]}x{manual_disp.shape[1]})', fontsize=9)\n", "for ax in axes:\n", " ax.axis('off')\n", "plt.colorbar(im, ax=axes, orientation='vertical', fraction=0.046, pad=0.04, label='disparity (px)')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "768b9d31", "metadata": {}, "source": [ "The border along the crop's edges can be ignored: the sliding window runs out of room to search there, an edge effect much like convolution's border handling (Lesson 10). Away from that border, and within the interior of the rectangles, the recovered disparity matches the ground truth almost perfectly, cleanly separating all three depth planes. The exception is the speckled band at each rectangle's edge: a window straddling a depth boundary mixes pixels from two different true disparities, so no single candidate disparity fits the whole window well, and the match becomes unreliable — a preview of the occlusion problem tackled next. Real photos are even less forgiving than this toy case: repetitive texture, flat regions, specular highlights, and shading all add genuine ambiguity on top of occlusion, which is why stereo matching is hard in practice." ] }, { "cell_type": "markdown", "id": "f80ab52d", "metadata": {}, "source": [ "## Catching bad matches: the left-right consistency check\n", "\n", "Occlusions and depth-boundary ambiguity produce wrong matches, like the speckle above. The fix: match both left-to-right and right-to-left, then keep only the pixels where the two results agree, discarding the rest as unreliable. Computing the right-to-left disparity needs no new code — flip both images left-right, swap which one plays \"left,\" run the exact same matcher, then flip the result back. A pixel fails the check whenever the two directions disagree by more than about a pixel. This is especially effective at catching occlusions: a pixel hidden behind a nearer surface in the other view has no true match there at all." ] }, { "cell_type": "code", "execution_count": null, "id": "f553f708", "metadata": {}, "outputs": [], "source": [ "block_size, max_disp = 7, 30\n", "crop_lr = np.s_[0:100, 0:100]\n", "left_c, right_c = left[crop_lr], right[crop_lr]\n", "\n", "disp_lr = block_match_stereo(left_c, right_c, block_size=block_size, max_disp=max_disp)\n", "disp_rl = block_match_stereo(right_c[:, ::-1], left_c[:, ::-1], block_size=block_size, max_disp=max_disp)[:, ::-1]\n", "\n", "h_c, w_c = disp_lr.shape\n", "consistent = np.zeros((h_c, w_c), dtype=bool)\n", "for y in range(h_c):\n", " for x in range(w_c):\n", " d = int(disp_lr[y, x])\n", " xr = x - d\n", " if 0 <= xr < w_c:\n", " consistent[y, x] = abs(disp_rl[y, xr] - d) <= 1\n", "\n", "# keep the full crop (no trimming) so the border-truncation effect is visible directly,\n", "# alongside the occlusion-driven inconsistency, in the figure below\n", "disp_checked = np.where(consistent, disp_lr, np.nan)\n", "\n", "fig, axes = plt.subplots(1, 3, figsize=(10, 3.5))\n", "for ax, d, title in zip(axes, [disp_lr, disp_rl, disp_checked],\n", " ['Left-to-right disparity', 'Right-to-left disparity', 'After consistency check']):\n", " im = ax.imshow(d, cmap=invalid_cmap, vmin=0, vmax=max_disp)\n", " ax.set_title(title, fontsize=9)\n", " ax.axis('off')\n", "plt.colorbar(im, ax=axes, orientation='vertical', fraction=0.046, pad=0.04, label='disparity (px)')\n", "plt.show()\n", "\n", "print(f'pixels marked inconsistent: {(~consistent).sum()} / {consistent.size}')" ] }, { "cell_type": "markdown", "id": "dc1a6c1c", "metadata": {}, "source": [ "Many pixels agree in both directions and survive the check. The ones that don't are typically caused by occlusion (the large gray regions next to vertical edges of the rectangles) or by the finite window size (small gray regions near the horizontal edges of the rectangles)." ] }, { "cell_type": "markdown", "id": "1da3f98f", "metadata": {}, "source": [ "## The same thing, at full resolution, with OpenCV\n", "\n", "`cv2.StereoBM` implements this same block-matching idea (with some efficiency and post-processing refinements) fast enough to run on the whole image." ] }, { "cell_type": "code", "execution_count": null, "id": "30fb738e", "metadata": {}, "outputs": [], "source": [ "stereo_bm = cv2.StereoBM_create(numDisparities=32, blockSize=9)\n", "disparity_bm = stereo_bm.compute(left, right).astype(np.float32) / 16.0 # fixed-point output, divide to get pixels\n", "\n", "print('recovered disparity vs. ground truth, by region:')\n", "print(f' background (true=5): {disparity_bm[100:110, 150:160].mean():.2f}')\n", "print(f' mid layer (true=15): {disparity_bm[40:50, 60:70].mean():.2f}')\n", "print(f' near layer (true=25): {disparity_bm[70:80, 90:110].mean():.2f}')\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(7, 3.5))\n", "im0 = axes[0].imshow(true_disparity, cmap='viridis', vmin=0, vmax=30)\n", "axes[0].set_title('Ground truth')\n", "im1 = axes[1].imshow(np.where(disparity_bm >= 0, disparity_bm, np.nan), cmap=invalid_cmap, vmin=0, vmax=30)\n", "axes[1].set_title('cv2.StereoBM output\\n(gray = invalid)')\n", "for ax in axes:\n", " ax.axis('off')\n", "plt.colorbar(im1, ax=axes, orientation='vertical', fraction=0.046, pad=0.04, label='disparity (px)')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "77d4eb0f", "metadata": {}, "source": [ "`StereoBM` marks a pixel invalid (returns $-1$) wherever it isn't confident in a unique match, e.g., too close to the image border to search the full disparity range. These pixels are shown in gray above." ] }, { "cell_type": "markdown", "id": "e210ffac", "metadata": {}, "source": [ "### Consistency check built into `disp12MaxDiff`\n", "\n", "The left-right consistency check is built into `cv2.StereoBM` — set `disp12MaxDiff` to a small positive number (typically 1) to enable it; the default of $-1$ leaves it off." ] }, { "cell_type": "code", "execution_count": null, "id": "05233606", "metadata": {}, "outputs": [], "source": [ "stereo_checked = cv2.StereoBM_create(numDisparities=32, blockSize=9)\n", "stereo_checked.setDisp12MaxDiff(1)\n", "disparity_checked = stereo_checked.compute(left, right).astype(np.float32) / 16.0\n", "\n", "newly_invalid = (disparity_bm >= 0) & (disparity_checked < 0)\n", "print(f'valid pixels without the check: {(disparity_bm >= 0).sum()}')\n", "print(f'valid pixels with the check: {(disparity_checked >= 0).sum()}')\n", "print(f'newly marked invalid: {newly_invalid.sum()}')\n", "\n", "plt.imshow(np.where(disparity_checked >= 0, disparity_checked, np.nan), cmap=invalid_cmap, vmin=0, vmax=30)\n", "plt.colorbar(label='disparity (px)')\n", "plt.title('cv2.StereoBM with disp12MaxDiff=1\\n(gray = invalid)')\n", "plt.axis('off')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "2c9c2784", "metadata": {}, "source": [ "## Window size: detail vs. noise tradeoff\n", "\n", "A larger matching window averages over more pixels, making the match more robust to noise but blurring across depth discontinuities (mixing pixels from two different true depths into one \"averaged\" disparity estimate near object edges). A smaller window preserves sharp depth boundaries but is more easily fooled by noise or repetitive texture." ] }, { "cell_type": "code", "execution_count": null, "id": "f2347d0e", "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(1, 3, figsize=(10, 3.5))\n", "for ax, block_size in zip(axes, [5, 15, 31]):\n", " stereo = cv2.StereoBM_create(numDisparities=32, blockSize=block_size)\n", " d = stereo.compute(left, right).astype(np.float32) / 16.0\n", " im = ax.imshow(np.where(d >= 0, d, np.nan), cmap=invalid_cmap, vmin=0, vmax=30)\n", " ax.set_title(f'blockSize={block_size}', fontsize=9)\n", " ax.axis('off')\n", "plt.colorbar(im, ax=axes, orientation='vertical', fraction=0.046, pad=0.04, label='disparity (px)')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "637a0fe5", "metadata": {}, "source": [ "## A real stereo pair: Tsukuba\n", "\n", "Everything so far has used synthetic images, precisely so the ground truth is known. For real photographs, the same three passes (left-to-right, right-to-left, and the consistency check between them) apply exactly as before, this time on an actual rectified stereo pair of photographs." ] }, { "cell_type": "markdown", "id": "9e45f41f", "metadata": {}, "source": [ "Flickering between the two frames makes the shift easy to see directly — watch how much farther the lamp moves than the bookshelves behind it, exactly the depth-dependent shift this lesson has been estimating all along.\n", "\n", "\"Flickering" ] }, { "cell_type": "markdown", "id": "519d065c", "metadata": {}, "source": [ "

Image source: Middlebury Stereo Datasets (University of Tsukuba)

" ] }, { "cell_type": "code", "execution_count": null, "id": "b3a4d7e2", "metadata": {}, "outputs": [], "source": [ "tsukuba_left = cv2.imread('../img/tsukuba_left.png', cv2.IMREAD_GRAYSCALE)\n", "tsukuba_right = cv2.imread('../img/tsukuba_right.png', cv2.IMREAD_GRAYSCALE)\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(7, 3.5))\n", "axes[0].imshow(tsukuba_left, cmap='gray')\n", "axes[0].set_title('Left photo', fontsize=9)\n", "axes[1].imshow(tsukuba_right, cmap='gray')\n", "axes[1].set_title('Right photo', fontsize=9)\n", "for ax in axes:\n", " ax.axis('off')\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "85212fd0", "metadata": {}, "outputs": [], "source": [ "num_disp = 64 # must be divisible by 16; large enough that the nearest surfaces don't clip the search range\n", "stereo = cv2.StereoBM_create(numDisparities=num_disp, blockSize=15)\n", "disp_lr = stereo.compute(tsukuba_left, tsukuba_right).astype(np.float32) / 16.0\n", "disp_rl = stereo.compute(tsukuba_right[:, ::-1], tsukuba_left[:, ::-1]).astype(np.float32) / 16.0\n", "disp_rl = disp_rl[:, ::-1]\n", "\n", "stereo_checked = cv2.StereoBM_create(numDisparities=num_disp, blockSize=15)\n", "stereo_checked.setDisp12MaxDiff(1)\n", "disp_checked = stereo_checked.compute(tsukuba_left, tsukuba_right).astype(np.float32) / 16.0\n", "\n", "changed = disp_lr != disp_checked\n", "print(f'pixels changed by the consistency check: {changed.sum()} / {changed.size} ({100 * changed.sum() / changed.size:.1f}%)')\n", "\n", "fig, axes = plt.subplots(1, 3, figsize=(10, 3.5))\n", "for ax, d, title in zip(axes, [disp_lr, disp_rl, disp_checked],\n", " ['Left-to-right disparity', 'Right-to-left disparity', 'After consistency check']):\n", " im = ax.imshow(np.where(d >= 0, d, np.nan), cmap=invalid_cmap, vmin=0, vmax=num_disp)\n", " ax.set_title(title, fontsize=9)\n", " ax.axis('off')\n", "plt.colorbar(im, ax=axes, orientation='vertical', fraction=0.046, pad=0.04, label='disparity (px)')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "29941d59", "metadata": {}, "source": [ "The lamp, the nearest object in the scene, stands out with the largest disparity; the bookshelves recede smoothly behind it. Only about 1% of pixels get flagged, so the third panel looks almost the same as the first at a glance — but that speckle is concentrated exactly along depth discontinuities throughout the scene: book spines, shelf edges, the head's silhouette, wherever a window straddles two different depths or one camera sees something the other doesn't." ] }, { "cell_type": "markdown", "id": "b0302e64", "metadata": {}, "source": [ "## From disparity to depth map\n", "\n", "Given a focal length and baseline (in whatever consistent units), $Z = fB/d$ converts a disparity map directly into a metric depth map." ] }, { "cell_type": "code", "execution_count": null, "id": "0466aa27", "metadata": {}, "outputs": [], "source": [ "focal_length_px = 500.0\n", "baseline_m = 0.1\n", "\n", "valid = disparity_bm > 0\n", "depth_m = np.full_like(disparity_bm, np.nan)\n", "depth_m[valid] = focal_length_px * baseline_m / disparity_bm[valid]\n", "\n", "plt.imshow(depth_m, cmap='viridis_r')\n", "plt.colorbar(label='estimated depth (m)')\n", "plt.title('Depth map (nearer = brighter)')\n", "plt.axis('off')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "a2130363", "metadata": {}, "outputs": [], "source": [ "step = 3 # subsample for a readable, fast-to-render point cloud\n", "ys, xs = np.mgrid[0:depth_m.shape[0]:step, 0:depth_m.shape[1]:step]\n", "zs = depth_m[::step, ::step]\n", "valid_pts = ~np.isnan(zs)\n", "\n", "fig = plt.figure(figsize=(6, 5))\n", "ax = fig.add_subplot(projection='3d')\n", "sc = ax.scatter(xs[valid_pts], ys[valid_pts], zs[valid_pts], c=zs[valid_pts], cmap='viridis_r', s=3)\n", "ax.set_xlabel('x (px)')\n", "ax.set_ylabel('y (px)')\n", "ax.set_zlabel('depth (m)')\n", "ax.invert_yaxis() # match image row convention (row 0 at top)\n", "ax.invert_zaxis() # nearer (smaller depth) points appear higher, like the disparity peaks above\n", "ax.view_init(elev=20, azim=-60)\n", "#fig.colorbar(sc, ax=ax, shrink=0.6, label='depth (m)')\n", "ax.set_title('The depth map, lifted into 3D')\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "ae597ee1", "metadata": {}, "source": [ "The same numbers, viewed as points floating in space instead of colors on a flat image, make the \"planes\" in this synthetic scene literal: three flat terraces at three different heights, with the holes where `StereoBM` had no confident match cut cleanly out of each one." ] }, { "cell_type": "markdown", "id": "d0f0011d", "metadata": {}, "source": [ "### Exercise\n", "\n", "1. Replace the random-texture `left` image with a flat gray region for the background (keeping the two rectangles textured). Run `cv2.StereoBM` again and describe what happens to the disparity estimate over the flat area — this is the aperture problem from Lesson 20, now in a stereo-matching context.\n", "2. Try `cv2.StereoSGBM_create` (semi-global block matching, which enforces smoothness across neighboring disparities rather than matching each pixel completely independently) in place of `StereoBM`. Compare the amount of speckle noise in flat/background regions between the two methods.\n", "3. Increase `baseline_m` in the depth conversion. How does the estimated depth map change, and why would a wider-baseline stereo rig give more precise depth estimates for distant objects (at the cost of a larger minimum-distance blind spot, due to disparity ranges and occlusion, near the cameras)?" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.x" } }, "nbformat": 4, "nbformat_minor": 5 }