{ "cells": [ { "cell_type": "markdown", "id": "2cd7620a", "metadata": {}, "source": "# Lesson 19: Feature Detection and Matching\n\nThis begins Part 2: classical 3D computer vision, where the algorithms increasingly work across *multiple* images of a scene — video frames, or several cameras — rather than a single one. The first problem that creates: two images of the same scene rarely put the same content at the same pixel location, since rotation, scale, and viewpoint all shift things around. This lesson builds toward feature-based matching, which survives exactly those changes: first **corner detectors** (Harris, Shi-Tomasi) that find distinctive, repeatable points, then **SIFT**, which adds scale-invariance and a descriptor that can be matched between two different-looking images of the same scene." }, { "cell_type": "code", "execution_count": null, "id": "d86f3619", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import cv2\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "a93774b9", "metadata": {}, "source": [ "## The structure tensor: what makes a good point to match?\n", "\n", "At each pixel, build the $2\\times2$ **structure tensor** from the local window of gradients (Lesson 12):\n", "\n", "$$M = \\sum_{\\text{window}} \\begin{bmatrix}I_x^2 & I_xI_y \\\\ I_xI_y & I_y^2\\end{bmatrix}$$\n", "\n", "This should look familiar: it's exactly the moment covariance matrix from Lesson 6, but built from gradients instead of pixel coordinates. Its eigenvalues $\\lambda_1 \\ge \\lambda_2$ describe the local intensity structure:\n", "\n", "- both small: **flat** region, gradients weak in every direction\n", "- one large, one small: **edge**, gradient strong perpendicular to the edge, weak along it\n", "- both large: **corner**, gradient strong in every direction — exactly the kind of point that can be precisely relocalized in another image" ] }, { "cell_type": "code", "execution_count": null, "id": "ad75f1cc", "metadata": {}, "outputs": [], "source": [ "img = np.zeros((100, 100), dtype=np.uint8)\n", "cv2.rectangle(img, (20, 20), (80, 80), 200, -1)\n", "cv2.rectangle(img, (40, 40), (60, 60), 100, -1)\n", "img_f = img.astype(np.float64)\n", "\n", "Ix = cv2.Sobel(img_f, cv2.CV_64F, 1, 0, ksize=3)\n", "Iy = cv2.Sobel(img_f, cv2.CV_64F, 0, 1, ksize=3)\n", "\n", "window = 5\n", "Sxx = cv2.boxFilter(Ix * Ix, -1, (window, window))\n", "Syy = cv2.boxFilter(Iy * Iy, -1, (window, window))\n", "Sxy = cv2.boxFilter(Ix * Iy, -1, (window, window))\n", "\n", "plt.imshow(img, cmap='gray')\n", "plt.title('Test image: corners, edges, and flat regions')\n", "plt.axis('off')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "cd745658", "metadata": {}, "source": [ "## Harris corner detector\n", "\n", "Directly computing eigenvalues everywhere is a bit expensive, so Harris and Stephens (1988) proposed a cheaper proxy using only $\\det(M)$ and $\\text{trace}(M)$ (which can be computed without ever finding eigenvalues):\n", "\n", "$$R = \\det(M) - k \\cdot \\text{trace}(M)^2, \\qquad k \\approx 0.04$$\n", "\n", "$R$ is large and positive at corners, negative at edges, and near zero on flat regions." ] }, { "cell_type": "code", "execution_count": null, "id": "5053eef9", "metadata": {}, "outputs": [], "source": [ "k = 0.04\n", "det_M = Sxx * Syy - Sxy**2\n", "trace_M = Sxx + Syy\n", "harris_response = det_M - k * trace_M**2\n", "\n", "harris_cv = cv2.cornerHarris(img_f.astype(np.float32), blockSize=window, ksize=3, k=k)\n", "print('correlation with cv2.cornerHarris:', np.corrcoef(harris_response.ravel(), harris_cv.ravel())[0, 1])\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(7, 3.5))\n", "axes[0].imshow(img, cmap='gray')\n", "axes[0].set_title('Image')\n", "im = axes[1].imshow(harris_response, cmap='coolwarm')\n", "axes[1].set_title('Harris response R\\n(red=corner, blue=edge)')\n", "for ax in axes:\n", " ax.axis('off')\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "f9cd5603", "metadata": {}, "source": [ "## Shi-Tomasi: \"Good Features to Track\"\n", "\n", "Shi and Tomasi (1994) argued for using the **minimum eigenvalue** of $M$ directly instead of Harris's determinant/trace proxy:\n", "\n", "$$R_{\\text{ST}} = \\min(\\lambda_1, \\lambda_2)$$\n", "\n", "A point is only a strong corner if *both* eigenvalues are large — the minimum being large is exactly that condition, and (unlike Harris's $R$) it has a direct, easily interpretable meaning: it's proportional to the worst-case tracking precision in any direction." ] }, { "cell_type": "code", "execution_count": null, "id": "bf4b74fb", "metadata": {}, "outputs": [], "source": [ "discriminant = np.sqrt(np.clip((trace_M / 2)**2 - det_M, 0, None))\n", "lambda_min = trace_M / 2 - discriminant\n", "\n", "eigval_cv = cv2.cornerMinEigenVal(img_f.astype(np.float32), blockSize=window, ksize=3)\n", "print('correlation with cv2.cornerMinEigenVal:', np.corrcoef(lambda_min.ravel(), eigval_cv.ravel())[0, 1])\n", "\n", "corners = cv2.goodFeaturesToTrack(img_f.astype(np.uint8), maxCorners=20, qualityLevel=0.1, minDistance=5)\n", "\n", "plt.imshow(img, cmap='gray')\n", "for pt in corners[:, 0]:\n", " plt.scatter(*pt, c='red', s=40, marker='+')\n", "plt.title('cv2.goodFeaturesToTrack (Shi-Tomasi)')\n", "plt.axis('off')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "2caf0562", "metadata": {}, "source": [ "## The scale problem\n", "\n", "Harris and Shi-Tomasi both operate at a single, fixed window size. Zoom the same corner in or out, and it may stop looking like a corner at that fixed scale (a sharp corner becomes a gentle curve when zoomed out enough) — these detectors are not scale-invariant. **SIFT** (Scale-Invariant Feature Transform) (Lowe, 1999/2004★) fixes this by searching for keypoints across an entire scale-space, not just one window size." ] }, { "cell_type": "markdown", "id": "6a8b4766", "metadata": {}, "source": [ "## SIFT: scale-space extrema and descriptors\n", "\n", "SIFT's pipeline, at a glance:\n", "\n", "1. Build a **Difference-of-Gaussians (DoG) scale-space** — exactly the DoG approximation to the Laplacian from Lesson 13, computed at many blur levels.\n", "2. Find keypoints as local extrema of DoG response across both space *and* scale — a point that's a maximum compared to its 26 neighbors (including 8 neighbors in the current scale *and* 9 neighbors in each scale level up and down).\n", "3. Assign each keypoint a dominant orientation from the local gradient histogram, so descriptors can be made rotation-invariant.\n", "4. Build a 128-dimensional descriptor from histograms of gradient orientations in a $4\\times4$ grid of subregions around the keypoint.\n", "\n", "The result: a keypoint with a position, a scale, an orientation, and a descriptor vector that's designed to be very similar between two images even under moderate rotation, scaling, and illumination change." ] }, { "cell_type": "code", "execution_count": null, "id": "c876173a", "metadata": {}, "outputs": [], "source": [ "textured = cv2.imread('../img/building.png', cv2.IMREAD_GRAYSCALE)\n", "textured = cv2.cvtColor(textured, cv2.COLOR_BGR2RGB)\n", "\n", "sift = cv2.SIFT_create()\n", "keypoints, descriptors = sift.detectAndCompute(textured, None)\n", "print(f'found {len(keypoints)} keypoints, each with a {descriptors.shape[1]}-dim descriptor')\n", "\n", "keypoint_vis = cv2.drawKeypoints(textured, keypoints, None,\n", " flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)\n", "\n", "fig, axes = plt.subplots(1, 2, figsize=(10, 5))\n", "axes[0].imshow(textured)\n", "axes[0].set_title('Original')\n", "axes[1].imshow(keypoint_vis)\n", "axes[1].set_title('SIFT keypoints (circle size = scale, line = orientation)')\n", "for ax in axes:\n", " ax.axis('off')\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "4d474d1d", "metadata": {}, "source": [ "

Photo by Dawson Lovell on Unsplash

" ] }, { "cell_type": "markdown", "id": "429094a0", "metadata": {}, "source": [ "## Matching under rotation and scale\n", "\n", "We rotate and shrink the image with a *known* transform, detect SIFT features independently in both, then match descriptors with a nearest-neighbor search. **Lowe's ratio test** keeps a match only if the best candidate is meaningfully closer than the *second*-best candidate — a simple, effective way to reject ambiguous matches." ] }, { "cell_type": "code", "execution_count": null, "id": "d3b08445", "metadata": {}, "outputs": [], "source": [ "M = cv2.getRotationMatrix2D((150, 150), angle=35, scale=0.7)\n", "rotated_scaled = cv2.warpAffine(textured, M, (300, 300))\n", "\n", "kp1, des1 = sift.detectAndCompute(textured, None)\n", "kp2, des2 = sift.detectAndCompute(rotated_scaled, 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)} (original), {len(kp2)} (rotated + scaled)')\n", "print(f'good matches after Lowe ratio test: {len(good_matches)} / {len(raw_matches)} candidate pairs')" ] }, { "cell_type": "markdown", "id": "bee9ac4a", "metadata": {}, "source": [ "### How many matches are actually correct?\n", "\n", "Because we know the exact transform `M` used to create the second image, we can directly check: does mapping each matched keypoint through `M` land near where it was actually matched?" ] }, { "cell_type": "code", "execution_count": null, "id": "12669abf", "metadata": {}, "outputs": [], "source": [ "correct = 0\n", "for m in good_matches:\n", " p1 = np.array(kp1[m.queryIdx].pt + (1,))\n", " p2 = np.array(kp2[m.trainIdx].pt)\n", " predicted = M @ p1\n", " if np.linalg.norm(predicted - p2) < 3:\n", " correct += 1\n", "\n", "print(f'geometrically correct matches: {correct} / {len(good_matches)} ({100 * correct / len(good_matches):.0f}%)')\n", "\n", "match_vis = cv2.drawMatches(textured, kp1, rotated_scaled, kp2, good_matches, None,\n", " flags=cv2.DrawMatchesFlags_NOT_DRAW_SINGLE_POINTS)\n", "plt.figure(figsize=(11, 5))\n", "plt.imshow(match_vis)\n", "plt.title(f'{len(good_matches)} matches across a 35-degree rotation + 0.7x scale change')\n", "plt.axis('off')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "c48c271c", "metadata": {}, "source": [ "Nearly all the ratio-test survivors are genuinely correct correspondences, despite the rotation and scale change — something a detector fixed to one scale and orientation simply couldn't recover. This is exactly what makes SIFT-style features the classic approach for tasks like image stitching, object recognition, and visual localization." ] }, { "cell_type": "markdown", "id": "05aa0df7", "metadata": {}, "source": [ "### Exercise\n", "\n", "1. Increase the rotation angle to 90 degrees and the scale to 0.3. Does the number of good matches and the geometric-correctness percentage hold up, or degrade? At what point does SIFT start to struggle?\n", "2. Try `cv2.ORB_create()` (a much faster, free binary-descriptor alternative to SIFT) with `cv2.NORM_HAMMING` in the `BFMatcher` instead of the default L2 norm. Compare the number of good matches and qualitatively compare speed using `%timeit` on `detectAndCompute`.\n", "3. Lower the ratio-test threshold from `0.75` to `0.6`. How do the number of good matches and the fraction that are geometrically correct both change? What does this tell you about the threshold's role in a precision/recall tradeoff?" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.x" } }, "nbformat": 4, "nbformat_minor": 5 }