{ "cells": [ { "cell_type": "markdown", "id": "e4d8059b", "metadata": {}, "source": "# Lesson 24: Image Stitching and Mosaicking\n\nThis lesson brings together feature detection and matching (Lesson 19), robust fitting (Lesson 22), and homography estimation (Lesson 23) into a complete pipeline that stitches two overlapping photos into a single seamless panorama." }, { "cell_type": "code", "execution_count": null, "id": "1536fa6b", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import cv2\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "533ed4ed", "metadata": {}, "source": [ "## The pipeline, at a glance\n", "\n", "1. **Detect and match features** between the two photos (Lesson 19: SIFT + ratio test).\n", "2. **Estimate a homography** relating one image's plane to the other's, robustly (Lesson 23's homography, Lesson 22's RANSAC).\n", "3. **Warp** one image into the other's coordinate frame (Lesson 9: `cv2.warpPerspective`).\n", "4. **Composite and blend** the two images onto one canvas, feathering across the overlap so the seam is invisible." ] }, { "cell_type": "markdown", "id": "91b0d1c4", "metadata": {}, "source": [ "## Two overlapping photos\n", "\n", "Two photos of the same brick building, taken from slightly different positions with substantial overlap." ] }, { "cell_type": "code", "execution_count": null, "id": "fa23d7f0", "metadata": {}, "outputs": [], "source": [ "view1 = cv2.imread('../img/clemson00.jpg')\n", "view2 = cv2.imread('../img/clemson01.jpg')\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(9, 3.5))\n", "axes[0].imshow(cv2.cvtColor(view1, cv2.COLOR_BGR2RGB))\n", "axes[0].set_title('View 1')\n", "axes[1].imshow(cv2.cvtColor(view2, cv2.COLOR_BGR2RGB))\n", "axes[1].set_title('View 2')\n", "for ax in axes:\n", " ax.axis('off')\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "45f2fee0", "metadata": {}, "source": [ "

Image source: Stan Birchfield

" ] }, { "cell_type": "markdown", "id": "61202b36", "metadata": {}, "source": [ "## Step 1: feature matching" ] }, { "cell_type": "code", "execution_count": null, "id": "3d5c8a18", "metadata": {}, "outputs": [], "source": [ "gray1 = cv2.cvtColor(view1, cv2.COLOR_BGR2GRAY)\n", "gray2 = cv2.cvtColor(view2, cv2.COLOR_BGR2GRAY)\n", "\n", "sift = cv2.SIFT_create()\n", "kp1, des1 = sift.detectAndCompute(gray1, None)\n", "kp2, des2 = sift.detectAndCompute(gray2, None)\n", "\n", "bf = cv2.BFMatcher()\n", "raw_matches = bf.knnMatch(des1, des2, k=2)\n", "good_matches = [m for m, n in raw_matches if m.distance < 0.75 * n.distance]\n", "\n", "print(f'keypoints: {len(kp1)} (view 1), {len(kp2)} (view 2)')\n", "print(f'good matches after ratio test: {len(good_matches)}')\n", "\n", "match_vis = cv2.drawMatches(view1, kp1, view2, kp2, good_matches[:60], None,\n", " flags=cv2.DrawMatchesFlags_NOT_DRAW_SINGLE_POINTS)\n", "plt.figure(figsize=(11, 4))\n", "plt.imshow(cv2.cvtColor(match_vis, cv2.COLOR_BGR2RGB))\n", "plt.title('Feature matches in the overlap region (first 60 shown)')\n", "plt.axis('off')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "ee1e85f3", "metadata": {}, "source": [ "## Step 2: robust homography\n", "\n", "We solve for the homography that maps points in view 2 into view 1's coordinate frame, so warping view 2 through it lands it in the right place on a shared canvas." ] }, { "cell_type": "code", "execution_count": null, "id": "2be36a3c", "metadata": {}, "outputs": [], "source": [ "pts1 = np.float32([kp1[m.queryIdx].pt for m in good_matches])\n", "pts2 = np.float32([kp2[m.trainIdx].pt for m in good_matches])\n", "\n", "H, inlier_mask = cv2.findHomography(pts2, pts1, cv2.RANSAC, 3.0)\n", "\n", "print(f'inliers: {int(inlier_mask.sum())} / {len(inlier_mask)}')\n", "print('recovered homography:')\n", "print(np.round(H, 4))" ] }, { "cell_type": "markdown", "id": "dee7e21b", "metadata": {}, "source": [ "## Steps 3-4: warp, composite, and blend\n", "\n", "First we figure out how big the output canvas needs to be, by mapping both images' corners through `H` (Lesson 23) and taking the bounding box — unlike the synthetic case, these photos aren't simple axis-aligned crops of a shared frame, so the canvas size and the offset of view 1 within it both have to be computed rather than assumed. We warp view 2 onto that canvas, then blend the overlap region with a simple **feather**: a linear alpha ramp from \"fully view 1\" to \"fully view 2\" across the overlap, the same weighted-sum blending as `cv2.addWeighted` in Lesson 2, just with a spatially-varying weight instead of a constant one." ] }, { "cell_type": "code", "execution_count": null, "id": "a00bb77a", "metadata": {}, "outputs": [], "source": [ "h1, w1 = view1.shape[:2]\n", "h2, w2 = view2.shape[:2]\n", "corners1 = np.float32([[0, 0], [w1, 0], [w1, h1], [0, h1]]).reshape(-1, 1, 2)\n", "corners2 = np.float32([[0, 0], [w2, 0], [w2, h2], [0, h2]]).reshape(-1, 1, 2)\n", "warped_corners2 = cv2.perspectiveTransform(corners2, H)\n", "all_corners = np.concatenate([corners1, warped_corners2], axis=0)\n", "\n", "x_min, y_min = np.floor(all_corners.min(axis=0).ravel()).astype(int)\n", "x_max, y_max = np.ceil(all_corners.max(axis=0).ravel()).astype(int)\n", "canvas_w, canvas_h = x_max - x_min, y_max - y_min\n", "\n", "translation = np.array([[1, 0, -x_min], [0, 1, -y_min], [0, 0, 1]], dtype=np.float64)\n", "canvas1 = np.zeros((canvas_h, canvas_w, 3), dtype=np.uint8)\n", "canvas1[-y_min:-y_min + h1, -x_min:-x_min + w1] = view1\n", "warped2 = cv2.warpPerspective(view2, translation @ H, (canvas_w, canvas_h))\n", "\n", "has1 = canvas1.sum(axis=2) > 0\n", "has2 = warped2.sum(axis=2) > 0\n", "overlap = has1 & has2\n", "\n", "overlap_cols = np.where(overlap.any(axis=0))[0]\n", "x_start, x_end = overlap_cols.min(), overlap_cols.max()\n", "ramp = np.clip((np.arange(canvas_w) - x_start) / (x_end - x_start + 1e-6), 0, 1)\n", "\n", "alpha = np.zeros((canvas_h, canvas_w), dtype=np.float32)\n", "alpha[overlap] = np.broadcast_to(ramp, (canvas_h, canvas_w))[overlap]\n", "\n", "stitched = canvas1.astype(np.float32) * (1 - alpha[..., None]) + warped2.astype(np.float32) * alpha[..., None]\n", "stitched[has1 & ~has2] = canvas1[has1 & ~has2] # regions covered only by view 1\n", "stitched[has2 & ~has1] = warped2[has2 & ~has1] # regions covered only by view 2\n", "stitched = stitched.astype(np.uint8)\n", "\n", "plt.figure(figsize=(10, 4))\n", "plt.imshow(cv2.cvtColor(stitched, cv2.COLOR_BGR2RGB))\n", "plt.title('Stitched panorama')\n", "plt.axis('off')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "609c526b", "metadata": {}, "source": [ "The roof (top, fully visible only in view 2) and the bush (bottom, fully visible only in view 1) both appear in the mosaic — confirming that the composite genuinely draws from both source images." ] }, { "cell_type": "markdown", "id": "8765d22f", "metadata": {}, "source": [ "### How good is the fit?\n", "\n", "There's no synthetic ground truth to compare against with a real photograph, but we don't need one: the RANSAC inliers themselves give a direct measure of geometric fit. For each inlier correspondence, project the point from view 2 through `H` and measure the distance to its matched point in view 1 — the **reprojection error**. A small, tightly clustered reprojection error means `H` explains the inlier geometry well." ] }, { "cell_type": "code", "execution_count": null, "id": "2bfb1ef0", "metadata": {}, "outputs": [], "source": [ "mask = inlier_mask.ravel().astype(bool)\n", "pts2_h = np.hstack([pts2[mask], np.ones((mask.sum(), 1))])\n", "proj = (H @ pts2_h.T).T\n", "proj = proj[:, :2] / proj[:, 2:3]\n", "reproj_error = np.linalg.norm(proj - pts1[mask], axis=1)\n", "\n", "print(f'mean reprojection error: {reproj_error.mean():.3f} px')\n", "print(f'max reprojection error: {reproj_error.max():.3f} px')\n", "\n", "plt.figure(figsize=(6, 3.5))\n", "plt.hist(reproj_error, bins=20)\n", "plt.xlabel('reprojection error (px)')\n", "plt.ylabel('inlier count')\n", "plt.title('Reprojection error of RANSAC inliers')\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "d590d554", "metadata": {}, "source": [ "Sub-pixel-scale reprojection error, well under the RANSAC inlier threshold of 3 px, confirms `H` is a good geometric fit — consistent with the clean seam visible in the stitched panorama above." ] }, { "cell_type": "markdown", "id": "68f47390", "metadata": {}, "source": [ "## In practice: `cv2.Stitcher`\n", "\n", "OpenCV bundles this entire pipeline (plus more robust blending, exposure compensation, and support for many images at once, arranged in any configuration) behind a single high-level call. Its default `PANORAMA` mode assumes the images come from a camera rotating about its optical center, and warps each image onto a sphere before compositing — great for wide panoramas, but it bows straight lines when the scene is dominated by flat, rectilinear structure (like a building facade) and the camera translated rather than purely rotated. `SCANS` mode instead composites with the planar homographies directly, matching the approach used above." ] }, { "cell_type": "code", "execution_count": null, "id": "7f8afda5", "metadata": {}, "outputs": [], "source": [ "stitcher = cv2.Stitcher_create(cv2.Stitcher_SCANS)\n", "status, panorama = stitcher.stitch([view1, view2])\n", "\n", "print('status:', 'OK' if status == cv2.Stitcher_OK else f'failed ({status})')\n", "if status == cv2.Stitcher_OK:\n", " plt.figure(figsize=(10, 4))\n", " plt.imshow(cv2.cvtColor(panorama, cv2.COLOR_BGR2RGB))\n", " plt.title('cv2.Stitcher result (SCANS mode)')\n", " plt.axis('off')\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "91673087", "metadata": {}, "source": [ "### Exercise\n", "\n", "1. Reduce the overlap between the two views (e.g. crop `view2` down to its rightmost quarter before matching). At what point does SIFT matching find too few good matches for `findHomography` to produce a reliable result?\n", "2. Replace the linear feather with a hard cutoff (no blending: just pick whichever image covers each pixel, splitting the overlap down the middle) and compare the visible seam quality to the feathered version.\n", "3. `H` here is a mild general homography, not a pure translation, because the two photos were taken from slightly different positions rather than a pure sideways pan. Print the ratio of `H`'s last row to `[0, 0, 1]` as a rough measure of how much perspective distortion it captures, and compare it to what you'd get by forcing an affine fit (`cv2.estimateAffinePartial2D`) instead — how well does the affine approximation stitch the images compared to the full homography?" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.x" } }, "nbformat": 4, "nbformat_minor": 5 }