Week 1 · Day 4 Lab — Edges, corners and descriptors¶

AIMS Cameroon · Introduction to Computer Vision

Learning goals

  1. Compute Sobel gradients and verify the worked example.
  2. Build Canny from parts: non-maximum suppression and hysteresis.
  3. Implement the Harris corner detector and test its invariances.
  4. Match SIFT and ORB features with Lowe's ratio test and measure how many matches are correct.
  5. Compute HOG and check its dimensionality.
In [ ]:
# If you are on Google Colab or Kaggle, everything below is pre-installed.
# Locally:  pip install numpy matplotlib scipy scikit-image scikit-learn opencv-python
import numpy as np
import matplotlib.pyplot as plt
import cv2, skimage, scipy
from skimage import data, color, util, exposure, filters, feature, measure, transform
from scipy import ndimage as ndi

plt.rcParams["figure.dpi"] = 100
rng = np.random.default_rng(0)

def show(*images, titles=None, cmap="gray", size=4, **kw):
    """Display images side by side."""
    n = len(images)
    fig, axes = plt.subplots(1, n, figsize=(size * n, size))
    axes = np.atleast_1d(axes)
    for i, (ax, im) in enumerate(zip(axes, images)):
        ax.imshow(im, cmap=cmap if np.ndim(im) == 2 else None, **kw)
        ax.axis("off")
        if titles:
            ax.set_title(titles[i])
    plt.tight_layout(); plt.show()

print("numpy", np.__version__, "| opencv", cv2.__version__, "| scikit-image", skimage.__version__)
In [ ]:
cam = data.camera().astype(float)
SOBEL_X = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], float)
SOBEL_Y = SOBEL_X.T

1. Gradients¶

✏️ Exercise 1.1¶

  1. Write gradients(img) returning Ix, Iy, mag, theta using ndi.correlate with the Sobel kernels (mode="reflect") and np.arctan2.
  2. Verify the slide example: the patch [[20,20,200]]*3 gives $G_x=720$, $G_y=0$.
In [ ]:
def gradients(img):
    # TODO
    raise NotImplementedError

patch = np.array([[20, 20, 200]] * 3, float)
# TODO: print Gx and Gy for the patch
In [ ]:
Ix, Iy, mag, theta = gradients(cam)
show(Ix, Iy, mag, theta, titles=["Ix", "Iy", "|∇I|", "θ"], size=3, cmap="gray")

2. Canny from parts¶

Steps: (1) Gaussian blur, (2) gradients, (3) non-maximum suppression, (4) double threshold, (5) hysteresis.

✏️ Exercise 2.1 — non-maximum suppression¶

Quantise the gradient direction (in degrees, modulo 180) to 0°, 45°, 90° or 135°. Keep a pixel only if its magnitude is ≥ both neighbours along the gradient direction:

direction neighbours (Δrow, Δcol)
0° (0, +1), (0, −1)
45° (+1, +1), (−1, −1)
90° (+1, 0), (−1, 0)
135° (+1, −1), (−1, +1)

(Rows point down in images, so with arctan2(Iy, Ix) these pairings are the correct ones. Check one by drawing it!)

In [ ]:
def nms(mag, theta):
    # TODO  hint: np.roll(mag, (-dr, -dc), axis=(0, 1))[r, c] == mag[r + dr, c + dc]
    raise NotImplementedError

✏️ Exercise 2.2 — double threshold + hysteresis¶

Strong: > hi. Weak: (lo, hi]. Keep weak pixels only if their 8-connected component (of strong ∪ weak pixels) contains at least one strong pixel. ndi.label(mask, structure=np.ones((3, 3))) does the connected components for you.

In [ ]:
def hysteresis(nms_mag, lo, hi):
    # TODO
    raise NotImplementedError

def my_canny(img, sigma=1.4, lo_pct=75, hi_pct=90):
    # TODO: blur -> gradients -> nms -> hysteresis (thresholds as percentiles of the magnitude)
    raise NotImplementedError
In [ ]:
mine = my_canny(cam)
ocv = cv2.Canny(cv2.GaussianBlur(data.camera(), (0, 0), 1.4), 50, 150) > 0
# tolerant agreement: an edge pixel "agrees" if the other map has an edge within 1 pixel
dil = lambda m: ndi.binary_dilation(m, np.ones((3, 3)))
print(f"my edges also in OpenCV (±1px): {(mine & dil(ocv)).sum() / mine.sum():.2f}")
print(f"OpenCV edges also in mine (±1px): {(ocv & dil(mine)).sum() / ocv.sum():.2f}")
show(mine, ocv, titles=["my_canny", "cv2.Canny"], size=5)

Questions. (a) Set lo_pct=0. What happens and why? (b) Set lo_pct=hi_pct. Which stage no longer matters? (c) Increase sigma to 3. What disappears first?

3. Harris corners¶

$M=\sum_W w\begin{bmatrix}I_x^2 & I_xI_y\\ I_xI_y& I_y^2\end{bmatrix}$, $R=\det M - k(\operatorname{tr}M)^2$.

✏️ Exercise 3.1¶

Implement harris(img, sigma=1.5, k=0.05) returning the response map $R$. Use a Gaussian window (cv2.GaussianBlur(..., (0, 0), sigma)) for $w$. Then detect corners with feature.peak_local_max(R, min_distance=6, threshold_rel=0.02).

In [ ]:
def harris(img, sigma=1.5, k=0.05):
    # TODO
    raise NotImplementedError
In [ ]:
board = data.checkerboard().astype(float)
R = harris(board)
pts = feature.peak_local_max(R, min_distance=6, threshold_rel=0.02)
print("corners found on the checkerboard:", len(pts))
plt.figure(figsize=(5, 5)); plt.imshow(board, cmap="gray"); plt.plot(pts[:, 1], pts[:, 0], "r+", ms=10); plt.axis("off"); plt.show()

✏️ Exercise 3.2 — invariances¶

First pad the board with a grey border (np.pad(board, 60, constant_values=127)) so that its outer corners are visible in every version. Then rotate it by 30° (transform.rotate(padded, 30, cval=127, preserve_range=True)), scale its brightness by 0.5, and downscale it by 4 (transform.rescale(..., preserve_range=True)). Count the corners each time (use the same threshold_rel). Which transformations is Harris (approximately) invariant to? Connect your answer to $M$.

In [ ]:
# TODO

4. SIFT and ORB matching¶

We create a second view of the astronaut with a known transform (rotation, scale, brightness), so we can check every match.

In [ ]:
ast = cv2.cvtColor(data.astronaut(), cv2.COLOR_RGB2GRAY)
h, w = ast.shape
M = cv2.getRotationMatrix2D((w / 2, h / 2), 30, 0.65)
view2 = cv2.warpAffine(ast, M, (w, h), borderValue=255)
view2 = np.clip(view2.astype(float) * 0.75 + 30, 0, 255).astype(np.uint8)
show(ast, view2, titles=["view 1", "view 2"])

✏️ Exercise 4.1¶

  1. Detect and describe with cv2.SIFT_create(); match with cv2.BFMatcher().knnMatch(d1, d2, k=2).
  2. Implement Lowe's ratio test (keep m if m.distance < 0.75 * n.distance).
  3. A match is correct if transforming the view-1 keypoint with M lands within 5 px of the view-2 keypoint. Report correct/total.
  4. Repeat with ORB (cv2.ORB_create(1000), cv2.BFMatcher(cv2.NORM_HAMMING)). Compare accuracy and speed.
In [ ]:
# TODO

5. HOG¶

✏️ Exercise 5.1¶

  1. Compute HOG on a 128×64 (rows×cols) pedestrian-sized window with 9 orientations, 8×8 cells, 2×2 blocks. Confirm the descriptor has 3780 dimensions, and derive this number on paper.
  2. Visualise HOG for the astronaut's face (visualize=True).
In [ ]:
# TODO

6. Stretch — a two-image panorama¶

Cut two overlapping crops from data.astronaut() (e.g. columns 0–320 and 180–512), detect and match SIFT features, estimate a homography with cv2.findHomography(src, dst, cv2.RANSAC, 3.0), and warp one crop onto the other with cv2.warpPerspective. How many RANSAC inliers do you get? What would break if the scene were not planar or the camera translated?

In [ ]:
# TODO