Lesson 24: Image Stitching and Mosaicking¶
This 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.
import numpy as np
import cv2
import matplotlib.pyplot as plt
The pipeline, at a glance¶
- Detect and match features between the two photos (Lesson 19: SIFT + ratio test).
- Estimate a homography relating one image's plane to the other's, robustly (Lesson 23's homography, Lesson 22's RANSAC).
- Warp one image into the other's coordinate frame (Lesson 9:
cv2.warpPerspective). - Composite and blend the two images onto one canvas, feathering across the overlap so the seam is invisible.
Two overlapping photos¶
Two photos of the same brick building, taken from slightly different positions with substantial overlap.
view1 = cv2.imread('../img/clemson00.jpg')
view2 = cv2.imread('../img/clemson01.jpg')
fig, axes = plt.subplots(1, 2, figsize=(9, 3.5))
axes[0].imshow(cv2.cvtColor(view1, cv2.COLOR_BGR2RGB))
axes[0].set_title('View 1')
axes[1].imshow(cv2.cvtColor(view2, cv2.COLOR_BGR2RGB))
axes[1].set_title('View 2')
for ax in axes:
ax.axis('off')
plt.tight_layout()
plt.show()
Image source: Stan Birchfield
Step 1: feature matching¶
gray1 = cv2.cvtColor(view1, cv2.COLOR_BGR2GRAY)
gray2 = cv2.cvtColor(view2, cv2.COLOR_BGR2GRAY)
sift = cv2.SIFT_create()
kp1, des1 = sift.detectAndCompute(gray1, None)
kp2, des2 = sift.detectAndCompute(gray2, None)
bf = cv2.BFMatcher()
raw_matches = bf.knnMatch(des1, des2, k=2)
good_matches = [m for m, n in raw_matches if m.distance < 0.75 * n.distance]
print(f'keypoints: {len(kp1)} (view 1), {len(kp2)} (view 2)')
print(f'good matches after ratio test: {len(good_matches)}')
match_vis = cv2.drawMatches(view1, kp1, view2, kp2, good_matches[:60], None,
flags=cv2.DrawMatchesFlags_NOT_DRAW_SINGLE_POINTS)
plt.figure(figsize=(11, 4))
plt.imshow(cv2.cvtColor(match_vis, cv2.COLOR_BGR2RGB))
plt.title('Feature matches in the overlap region (first 60 shown)')
plt.axis('off')
plt.show()
Step 2: robust homography¶
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.
pts1 = np.float32([kp1[m.queryIdx].pt for m in good_matches])
pts2 = np.float32([kp2[m.trainIdx].pt for m in good_matches])
H, inlier_mask = cv2.findHomography(pts2, pts1, cv2.RANSAC, 3.0)
print(f'inliers: {int(inlier_mask.sum())} / {len(inlier_mask)}')
print('recovered homography:')
print(np.round(H, 4))
Steps 3-4: warp, composite, and blend¶
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.
h1, w1 = view1.shape[:2]
h2, w2 = view2.shape[:2]
corners1 = np.float32([[0, 0], [w1, 0], [w1, h1], [0, h1]]).reshape(-1, 1, 2)
corners2 = np.float32([[0, 0], [w2, 0], [w2, h2], [0, h2]]).reshape(-1, 1, 2)
warped_corners2 = cv2.perspectiveTransform(corners2, H)
all_corners = np.concatenate([corners1, warped_corners2], axis=0)
x_min, y_min = np.floor(all_corners.min(axis=0).ravel()).astype(int)
x_max, y_max = np.ceil(all_corners.max(axis=0).ravel()).astype(int)
canvas_w, canvas_h = x_max - x_min, y_max - y_min
translation = np.array([[1, 0, -x_min], [0, 1, -y_min], [0, 0, 1]], dtype=np.float64)
canvas1 = np.zeros((canvas_h, canvas_w, 3), dtype=np.uint8)
canvas1[-y_min:-y_min + h1, -x_min:-x_min + w1] = view1
warped2 = cv2.warpPerspective(view2, translation @ H, (canvas_w, canvas_h))
has1 = canvas1.sum(axis=2) > 0
has2 = warped2.sum(axis=2) > 0
overlap = has1 & has2
overlap_cols = np.where(overlap.any(axis=0))[0]
x_start, x_end = overlap_cols.min(), overlap_cols.max()
ramp = np.clip((np.arange(canvas_w) - x_start) / (x_end - x_start + 1e-6), 0, 1)
alpha = np.zeros((canvas_h, canvas_w), dtype=np.float32)
alpha[overlap] = np.broadcast_to(ramp, (canvas_h, canvas_w))[overlap]
stitched = canvas1.astype(np.float32) * (1 - alpha[..., None]) + warped2.astype(np.float32) * alpha[..., None]
stitched[has1 & ~has2] = canvas1[has1 & ~has2] # regions covered only by view 1
stitched[has2 & ~has1] = warped2[has2 & ~has1] # regions covered only by view 2
stitched = stitched.astype(np.uint8)
plt.figure(figsize=(10, 4))
plt.imshow(cv2.cvtColor(stitched, cv2.COLOR_BGR2RGB))
plt.title('Stitched panorama')
plt.axis('off')
plt.show()
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.
How good is the fit?¶
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.
mask = inlier_mask.ravel().astype(bool)
pts2_h = np.hstack([pts2[mask], np.ones((mask.sum(), 1))])
proj = (H @ pts2_h.T).T
proj = proj[:, :2] / proj[:, 2:3]
reproj_error = np.linalg.norm(proj - pts1[mask], axis=1)
print(f'mean reprojection error: {reproj_error.mean():.3f} px')
print(f'max reprojection error: {reproj_error.max():.3f} px')
plt.figure(figsize=(6, 3.5))
plt.hist(reproj_error, bins=20)
plt.xlabel('reprojection error (px)')
plt.ylabel('inlier count')
plt.title('Reprojection error of RANSAC inliers')
plt.tight_layout()
plt.show()
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.
In practice: cv2.Stitcher¶
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.
stitcher = cv2.Stitcher_create(cv2.Stitcher_SCANS)
status, panorama = stitcher.stitch([view1, view2])
print('status:', 'OK' if status == cv2.Stitcher_OK else f'failed ({status})')
if status == cv2.Stitcher_OK:
plt.figure(figsize=(10, 4))
plt.imshow(cv2.cvtColor(panorama, cv2.COLOR_BGR2RGB))
plt.title('cv2.Stitcher result (SCANS mode)')
plt.axis('off')
plt.show()
Exercise¶
- Reduce the overlap between the two views (e.g. crop
view2down to its rightmost quarter before matching). At what point does SIFT matching find too few good matches forfindHomographyto produce a reliable result? - 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.
Hhere 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 ofH'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?