CSCI 3240U: Computer Vision I

Lab 7: Interest Points and Local Descriptors

Faisal Z. Qureshi
Faculty of Science, Ontario Tech University
Oshawa ON Canada
http://vclab.science.ontariotechu.ca

Fall 2026

Copyright information

© Faisal Qureshi

License

This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.

About this notebook

This notebook is the starter for Part 1 of Lab 7, corner detection. Read the handout first. The sections below follow its numbered tasks.

Section here	Handout
Task 1: Harris corner detector	Part 1, Task 1
Non-maximum suppression	Deliverables (no numbered task asks for it, but the Deliverables do)
Task 2: Shi-Tomasi detector	Part 1, Task 2
Task 3: compare the two detectors	Part 1, Task 3
Task 4: rotation repeatability	Part 1, Task 4
Optional extra: img_match	Part 4, optional and not marked

The other parts of the lab live elsewhere. Part 2 (optical flow and tracking) has no starter notebook; write it yourself against ../lab1-setup/traffic-short.mp4. Part 3 (the SIFT-like descriptor) has its own starter, sift-like-descriptor.ipynb.

Cells marked TODO are yours to write. Everything else runs as given and is there to save you plumbing.

Submit a single notebook, executed top to bottom, via Canvas.

Setup

import os

import cv2
import numpy as np
import matplotlib
import matplotlib.pyplot as plt

%matplotlib inline

print('OpenCV', cv2.__version__)
print('NumPy', np.__version__)
print('Matplotlib', matplotlib.__version__)

Helper code

Nothing below is graded. These helpers load images, draw points, rotate an image by a known angle, move point coordinates between the original and the rotated frame, and plot a repeatability curve. Read them before you use them.

Two conventions matter and they are easy to mix up.

• OpenCV point coordinates are (x, y), that is (column, row).
• NumPy array indices are [row, column].

Every point array in this notebook is (N, 2) holding (x, y), in float.

DATA_FOLDER = '.'


def load_rgb(filename):
    """Read an image and return it as RGB, ready for plt.imshow.

    cv2.imread gives you BGR. Passing that straight to plt.imshow swaps red and
    blue, which is the classic first bug of every OpenCV lab.
    """
    path = os.path.join(DATA_FOLDER, filename)
    bgr = cv2.imread(path)
    if bgr is None:
        raise FileNotFoundError(path)
    return cv2.cvtColor(bgr, cv2.COLOR_BGR2RGB)


def load_gray(filename):
    """Read an image as a single channel uint8 array."""
    path = os.path.join(DATA_FOLDER, filename)
    gray = cv2.imread(path, cv2.IMREAD_GRAYSCALE)
    if gray is None:
        raise FileNotFoundError(path)
    return gray


def show(img, title=None, size=(10, 10), cmap=None):
    """Display one image.  Pass an RGB array, not a BGR one."""
    plt.figure(figsize=size)
    if title:
        plt.title(title)
    if img.ndim == 2 and cmap is None:
        cmap = 'gray'
    plt.imshow(img, cmap=cmap)
    plt.axis('off')


def draw_points(img, pts, colour=(255, 0, 0), radius=4, thickness=-1):
    """Draw (x, y) points on a copy of an RGB image and return the copy."""
    out = img.copy()
    if pts is None:
        return out
    for x, y in np.asarray(pts, dtype=float).reshape(-1, 2):
        cv2.circle(out, (int(round(x)), int(round(y))), radius, colour, thickness)
    return out

def rotate_image(img, angle_deg):
    """Rotate about the image centre onto a canvas large enough to hold it all.

    Returns the rotated image and the 2x3 affine matrix M that maps a point
    in the original image to its location in the rotated image.
    """
    h, w = img.shape[:2]
    centre = (w / 2.0, h / 2.0)
    M = cv2.getRotationMatrix2D(centre, angle_deg, 1.0)

    cos, sin = abs(M[0, 0]), abs(M[0, 1])
    new_w = int(round(h * sin + w * cos))
    new_h = int(round(h * cos + w * sin))

    # Shift so the rotated content lands inside the new canvas.
    M[0, 2] += new_w / 2.0 - centre[0]
    M[1, 2] += new_h / 2.0 - centre[1]

    rotated = cv2.warpAffine(img, M, (new_w, new_h), flags=cv2.INTER_LINEAR)
    return rotated, M


def transform_points(pts, M):
    """Apply a 2x3 affine matrix to an (N, 2) array of (x, y) points."""
    pts = np.asarray(pts, dtype=np.float64).reshape(-1, 2)
    if len(pts) == 0:
        return pts
    homogeneous = np.hstack([pts, np.ones((len(pts), 1))])
    return homogeneous @ M.T


def inverse_affine(M):
    """Invert a 2x3 affine matrix.  Use this to map rotated-frame points back."""
    return cv2.invertAffineTransform(M)


def inside(pts, shape, margin=0):
    """Boolean mask: which (x, y) points fall inside an image of this shape."""
    pts = np.asarray(pts, dtype=float).reshape(-1, 2)
    h, w = shape[:2]
    return ((pts[:, 0] >= margin) & (pts[:, 0] < w - margin) &
            (pts[:, 1] >= margin) & (pts[:, 1] < h - margin))


def plot_repeatability(angles, curves, tol):
    """Plot repeatability against rotation angle.

    curves is a dict mapping a label to a list of rates, one per angle.
    """
    plt.figure(figsize=(7, 4))
    for label, rates in curves.items():
        plt.plot(angles, rates, marker='o', label=label)
    plt.xlabel('rotation (degrees)')
    plt.ylabel('repeatability')
    plt.title('Repeatability within %g px' % tol)
    plt.ylim(0.0, 1.05)
    plt.grid(alpha=0.3)
    plt.legend()

Test image

cn-tower-1.jpg is a reasonable default. cb.png is a synthetic checkerboard and is useful while you are debugging, because you know where its corners are.

img_file = 'cn-tower-1.jpg'
# img_file = 'cb.png'

img = load_rgb(img_file)
gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)

print('image shape', img.shape)
show(img, 'Test image')

Task 1: Harris corner detector

Both detectors in this lab are built from the structure tensor

$$ M = \sum_{(x,y) \in W} w(x,y) \left[ \begin{array}{cc} I_x^2 & I_x I_y \\ I_x I_y & I_y^2 \end{array} \right], $$

and they differ only in how they turn its eigenvalues into a score. Harris uses $R = \det(M) - k\,\mathrm{tr}(M)^2$ with $k$ around 0.04 to 0.06.

The cell below runs OpenCV's cornerHarris and shows the raw response. Use it as a reference. You may also build $M$ yourself from the $I_x$ and $I_y$ machinery of Lab 4; the two responses should agree up to a scale factor.

src = np.float32(gray)

blocksize = 2   # neighbourhood used to accumulate the structure tensor
ksize = 3       # Sobel aperture used for the derivatives
k = 0.04        # Harris parameter in R = det(M) - k (trace M)^2

harris_response = cv2.cornerHarris(src, blocksize, ksize, k)

# A crude threshold on the raw response. Note the blobs: each real corner fires
# over a small neighbourhood, so one corner produces many detections. That is
# the problem non-maximum suppression solves.
crude = img.copy()
crude[harris_response > 0.01 * harris_response.max()] = [255, 0, 0]

plt.figure(figsize=(20, 15))
plt.subplot(2, 1, 1)
plt.title('Harris response R')
plt.imshow(harris_response, cmap='gray')
plt.axis('off')
plt.subplot(2, 1, 2)
plt.title('Thresholded response, no suppression')
plt.imshow(crude)
plt.axis('off')

print('pixels above threshold:', int((harris_response > 0.01 * harris_response.max()).sum()))

TODO (Task 1)

Turn the response map into a list of interest points.

Write harris_keypoints. It must return an (N, 2) float array of (x, y) locations, sorted by response so that max_points keeps the strongest ones. Threshold the response relative to its maximum, then suppress non-maxima with the nms function you write in the next section.

def harris_keypoints(gray_img, max_points=300, rel_threshold=0.01,
                     blocksize=2, ksize=3, k=0.04, nms_window=7):
    """Return an (N, 2) array of Harris interest points as (x, y)."""
    # TODO
    #   1. compute the response (cv2.cornerHarris, or your own structure tensor)
    #   2. threshold it at rel_threshold * response.max()
    #   3. apply non-maximum suppression over an nms_window neighbourhood
    #   4. sort the survivors by response, strongest first, and keep max_points
    raise NotImplementedError('Task 1: harris_keypoints')

Non-maximum suppression

The Deliverables require your detections to be non-maximum suppressed. No numbered task says so in as many words, so it is easy to miss. Without it you report a blob of detections per corner, and every count you make afterwards, including the Task 4 repeatability, is meaningless.

The idea is one line: keep a pixel only if it is the largest response in its own neighbourhood. A max filter over the response makes this easy. cv2.dilate with a rectangular structuring element is a max filter, so response == cv2.dilate(response, kernel) marks the local maxima. Combine that with your threshold.

Note that cv2.dilate used on its own, as many online tutorials do, does the opposite of suppression: it spreads each response outwards so the marks are easier to see. Do not confuse the two uses.

def nms(response, threshold, window=7):
    """Non-maximum suppression on a response map.

    Keep a pixel if its response exceeds `threshold` and is the maximum inside
    a window x window neighbourhood centred on it.

    Returns an (N, 2) float array of (x, y) locations.
    """
    # TODO
    raise NotImplementedError('non-maximum suppression')

# Check your suppression once nms is written. Both counts print, so you can see
# how much the suppression removed. Expect a large reduction.
thr = 0.01 * harris_response.max()
raw_count = int((harris_response > thr).sum())
pts_nms = nms(harris_response, thr, window=7)

print('above threshold      :', raw_count)
print('after suppression    :', len(pts_nms))

plt.figure(figsize=(20, 8))
plt.subplot(1, 2, 1)
plt.title('No suppression')
plt.imshow(crude)
plt.axis('off')
plt.subplot(1, 2, 2)
plt.title('After suppression')
plt.imshow(draw_points(img, pts_nms, colour=(255, 0, 0), radius=4))
plt.axis('off')

Task 2: Shi-Tomasi detector

Shi and Tomasi score a window by $R = \min(\lambda_1, \lambda_2)$, the smaller eigenvalue of the same $M$. Part 2 of this lab explains why that particular score earns the title Good Features to Track.

cv2.goodFeaturesToTrack implements it. Its minDistance argument performs the suppression for you, so its output is already sparse. Keep that in mind when you compare it against Harris: compare like with like, by giving both detectors the same suppression window and the same point budget.

corners = cv2.goodFeaturesToTrack(gray,
                                  maxCorners=25,
                                  qualityLevel=0.01,
                                  minDistance=10)

# np.int0 was removed in NumPy 2.0. Use np.intp, or just keep the floats and
# round when you draw. Here we keep floats, which is what the helpers expect.
pts_shi_demo = corners.reshape(-1, 2).astype(np.float64)

print('Shi-Tomasi points:', len(pts_shi_demo))
show(draw_points(img, pts_shi_demo, colour=(0, 255, 0), radius=5),
     'Shi-Tomasi corner detector')

TODO (Task 2)

Write shi_tomasi_keypoints with the same signature style as harris_keypoints, so that Task 3 can call both the same way and compare them fairly. Return an (N, 2) float array of (x, y) locations.

You may call cv2.goodFeaturesToTrack, or compute $\min(\lambda_1,\lambda_2)$ yourself with cv2.cornerMinEigenVal and reuse your nms. If you use goodFeaturesToTrack, set minDistance to match your Harris suppression window, otherwise the comparison in Task 3 measures your parameters rather than the detectors.

def shi_tomasi_keypoints(gray_img, max_points=300, rel_threshold=0.01,
                         nms_window=7):
    """Return an (N, 2) array of Shi-Tomasi interest points as (x, y)."""
    # TODO
    raise NotImplementedError('Task 2: shi_tomasi_keypoints')

Task 3: compare the two detectors

Run both detectors on the same image with matched settings, draw both point sets on one figure, and say what differs.

MAX_POINTS = 300
NMS_WINDOW = 7
REL_THRESHOLD = 0.01

pts_harris = harris_keypoints(gray, max_points=MAX_POINTS,
                              rel_threshold=REL_THRESHOLD,
                              nms_window=NMS_WINDOW)
pts_shi = shi_tomasi_keypoints(gray, max_points=MAX_POINTS,
                               rel_threshold=REL_THRESHOLD,
                               nms_window=NMS_WINDOW)

print('Harris     :', len(pts_harris), 'points')
print('Shi-Tomasi :', len(pts_shi), 'points')

overlay = draw_points(img, pts_harris, colour=(255, 0, 0), radius=6, thickness=1)
overlay = draw_points(overlay, pts_shi, colour=(0, 255, 0), radius=3, thickness=-1)
show(overlay, 'Harris (red circles) and Shi-Tomasi (green dots)', size=(14, 14))

TODO (Task 3)

"They look about the same" is a fair first impression, and it is not an answer. Quantify it.

Write a function that, given the two point sets and a tolerance in pixels, reports what fraction of Harris points have a Shi-Tomasi point within that tolerance, and the fraction the other way round. Then write a short paragraph in the markdown cell below: where do the two agree, where do they disagree, and what kind of image structure sits at the points only one of them found?

def overlap_fraction(pts_a, pts_b, tol=3.0):
    """Fraction of points in pts_a with a point of pts_b within tol pixels."""
    # TODO
    raise NotImplementedError('Task 3: overlap_fraction')

Your Task 3 write-up goes here.

Task 4: rotation repeatability

Rotate the image by a known angle, detect again, and ask whether the same physical points come back.

The measurement has to be made on locations, not on counts. A detector that finds 300 points before the rotation and 300 after has told you nothing: they could be 300 completely different points. So:

1. detect in the original image, giving pts_ref;
2. rotate the image by $\theta$ with rotate_image, which also hands you the affine matrix $M$;
3. detect in the rotated image, giving pts_rot;
4. map pts_rot back into the original frame with the inverse of $M$;
5. count how many of pts_ref have a mapped-back point within a few pixels.

Step 4 is pure geometry and is provided. Step 5 is the measurement and is yours.

One detail decides whether your number means anything. rotate_image enlarges the canvas, so the rotated image carries blank corners and a long straight border that was never in the scene. Detections there are artefacts. Keep only the mapped-back points that land inside the original image, and be ready to discard a border margin as well. inside() does both.

ANGLE = 30.0

rot_img, M = rotate_image(img, ANGLE)
rot_gray = cv2.cvtColor(rot_img, cv2.COLOR_RGB2GRAY)
M_inv = inverse_affine(M)

print('original shape', img.shape[:2], ' rotated shape', rot_img.shape[:2])

# Sanity check on the geometry: send a point out and bring it back.
probe = np.array([[100.0, 200.0], [400.0, 50.0]])
round_trip = transform_points(transform_points(probe, M), M_inv)
print('round trip error (px):', np.abs(round_trip - probe).max())

show(rot_img, 'Rotated by %g degrees' % ANGLE, size=(12, 12))

TODO (Task 4)

Write repeatability. Given the reference points and the mapped-back points, return the fraction of reference points that were recovered.

Decide, and state in your write-up, how you handle these:

• a mapped-back point should be used at most once, otherwise one lucky detection can claim several reference points;
• a reference point sitting on the image border loses part of its neighbourhood once the image is rotated, so decide whether to exclude a margin;
• what tolerance you use, and why. A few pixels is the usual choice.

The driver cell below already drops mapped-back points that fall outside the original image, so repeatability receives two point sets in the same frame.

def repeatability(pts_ref, pts_mapped, tol=3.0):
    """Fraction of pts_ref matched by a distinct point of pts_mapped within tol."""
    # TODO
    raise NotImplementedError('Task 4: repeatability')

# Driver. This plumbing is given; the number it reports comes from your
# repeatability function, so it is only as good as that.
TOL = 3.0
angles = [0, 10, 20, 30, 45, 60, 90]
curves = {'Harris': [], 'Shi-Tomasi': []}

detectors = {'Harris': harris_keypoints, 'Shi-Tomasi': shi_tomasi_keypoints}

for name, detect in detectors.items():
    ref = detect(gray, max_points=MAX_POINTS, rel_threshold=REL_THRESHOLD,
                 nms_window=NMS_WINDOW)
    for angle in angles:
        rotated, Ma = rotate_image(img, angle)
        rotated_gray = cv2.cvtColor(rotated, cv2.COLOR_RGB2GRAY)
        found = detect(rotated_gray, max_points=MAX_POINTS,
                       rel_threshold=REL_THRESHOLD, nms_window=NMS_WINDOW)
        back = transform_points(found, inverse_affine(Ma))
        back = back[inside(back, img.shape)]   # drop points off the original canvas
        curves[name].append(repeatability(ref, back, tol=TOL))
    print(name, ['%.2f' % r for r in curves[name]])

plot_repeatability(angles, curves, TOL)

# Look at the 30 degree case. Reference points in red, points detected in the
# rotated image and mapped back in green. Where the two coincide the detector
# was repeatable.
ref_h = harris_keypoints(gray, max_points=MAX_POINTS,
                         rel_threshold=REL_THRESHOLD, nms_window=NMS_WINDOW)
rot_h = harris_keypoints(rot_gray, max_points=MAX_POINTS,
                         rel_threshold=REL_THRESHOLD, nms_window=NMS_WINDOW)
back_h = transform_points(rot_h, M_inv)
back_h = back_h[inside(back_h, img.shape)]

vis = draw_points(img, ref_h, colour=(255, 0, 0), radius=6, thickness=1)
vis = draw_points(vis, back_h, colour=(0, 255, 0), radius=3, thickness=-1)
show(vis, 'Original detections (red) and rotated detections mapped back (green)',
     size=(14, 14))

Your Task 4 write-up goes here. Report the repeatability at 30 degrees for both detectors, with the tolerance you used. Then account for the points that were lost. Interpolation, the discrete detection window and the image boundary all cost you something. If your number is close to 100 percent, check that you compared locations rather than counts.

Part 4 of the handout: deciding whether two images match

This part is optional and is not marked. Everything above is. Come back to it if you have time, or after Lab 8, which does descriptor matching properly inside RANSAC.

SIFT keypoints

The code below finds locations suitable for SIFT, each with a scale and an orientation.

SIFT moved into the main OpenCV module in version 4.4. Call cv2.SIFT_create(). The old cv2.xfeatures2d.SIFT_create() does not exist in opencv-python and raises an AttributeError.

sift = cv2.SIFT_create()
kp = sift.detect(gray, None)

print('SIFT keypoints:', len(kp))

# Draw each keypoint with its scale and orientation.
sift_vis = cv2.drawKeypoints(gray, kp, None,
                             flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
show(sift_vis, 'SIFT key points')

TODO (optional, not marked)

Use local feature descriptors from two images to compute a match score for the pair. Here is a recipe. You are encouraged not to use the built-in matcher. Experiment with different matching schemes instead.

1. Compute $n$ SIFT descriptors from image 1.
2. Compute $m$ SIFT descriptors from image 2.
3. Each descriptor is a 128-dimensional vector, so pick a distance between two vectors, say Euclidean, and a threshold that decides whether two vectors match. Using that distance, count the matched vectors between the two sets.
4. Report the matched count as a fraction of the smaller of $n$ and $m$.

Complete the following function, which returns a number between 0.0 and 1.0.

def img_match(filename1, filename2):

    # TO DO

    return 0.0

We will use this function as follows.

filename1 = 'cn-tower-1.jpg'
filename2 = 'cn-tower-2.jpg'

print ('Match score is', img_match(filename1, filename2))

Submission

Submit a single Jupyter notebook via Canvas. Run it top to bottom before you submit, so every figure and every number is visible in the file. Code that has not been executed cannot be marked.