Corners & Features: Harris & ORB

By the end of this lesson you will understand how a computer finds the same point in two different photos. You will code the Harris corner detector, build a tiny binary descriptor, compare descriptors with Hamming distance, and use OpenCV's ORB to match features between a photo and a rotated copy. This is how panorama apps stitch photos together.

Why corners?

Suppose you take two photos of the Hawa Mahal in Jaipur, the second one a step to the left. To stitch them into a panorama, the computer must find the same points in both photos. Which points are easy to find again?

• A point in a flat area (blue sky) looks like every other sky point. Useless.
• A point on an edge (the roof line) can slide along the edge and still look the same. Its position is uncertain along the edge.
• A corner (the tip of a window frame) looks different if you move it in any direction. It can be located exactly.

That is why corners, also called keypoints or interest points, are the starting point of feature matching. Think of giving directions in a crowded market: "near the wall" is vague, but "at the corner where the sweet shop meets the chemist" is exact.

The Harris corner detector

Harris measures how much the image changes when a small window shifts. It uses the gradients gx and gy from Edge Detection, collected in a small neighbourhood around each pixel:

• Sxx = sum of gx², Syy = sum of gy², Sxy = sum of gx × gy (summed with a Gaussian window).

These three numbers form a 2 by 2 matrix called the structure tensor. Its two eigenvalues (two numbers that describe how strongly the brightness changes in the two main directions) tell the story:

Computing eigenvalues for every pixel is slow, so Harris uses a shortcut score that needs no eigenvalues:

R = (Sxx × Syy - Sxy²) - k × (Sxx + Syy)², with k about 0.04 to 0.06.

R is large and positive at corners, negative on edges and near zero in flat areas.

Step 1: compute gradients. Step 2: build the three sums with a Gaussian window. Step 3: compute R. Step 4: keep the local maxima above a threshold.

import numpy as np
import matplotlib.pyplot as plt
from numpy.lib.stride_tricks import sliding_window_view

def conv2d(img, kernel):
    p = np.pad(img.astype(np.float64), kernel.shape[0] // 2, mode="edge")
    return np.einsum("ijkl,kl->ij", sliding_window_view(p, kernel.shape), kernel)

def gaussian_kernel(size, sigma):
    h = size // 2
    y, x = np.mgrid[-h:h + 1, -h:h + 1]
    k = np.exp(-(x ** 2 + y ** 2) / (2 * sigma ** 2))
    return k / k.sum()

img = np.full((100, 140), 40.0)
img[20:60, 20:70] = 200                      # a rectangle: 4 corners
img[50:90, 90:130] = 140                     # a square: 4 more corners

# Step 1
sx = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
gx, gy = conv2d(img, sx), conv2d(img, sx.T)
# Step 2
g = gaussian_kernel(5, 1.0)
Sxx, Syy, Sxy = conv2d(gx * gx, g), conv2d(gy * gy, g), conv2d(gx * gy, g)
# Step 3
k = 0.05
R = (Sxx * Syy - Sxy ** 2) - k * (Sxx + Syy) ** 2
# Step 4: local maximum in a 7x7 window and above 1% of the strongest response
p = np.pad(R, 3, constant_values=-np.inf)
local_max = R == sliding_window_view(p, (7, 7)).max(axis=(2, 3))
corners = np.argwhere(local_max & (R > 0.01 * R.max()))
print("corners found:", len(corners))
print(corners)
print(f"R at a corner: {R[20, 20]:.3g}   on an edge: {R[20, 45]:.3g}   in a flat area: {R[80, 30]:.3g}")

plt.imshow(img, cmap="gray")
plt.scatter(corners[:, 1], corners[:, 0], c="red", s=30)
plt.title("Harris corners")
plt.show()

corners found: 8
[[ 20  20]
 [ 20  69]
 [ 50  90]
 [ 50 129]
 [ 59  20]
 [ 59  69]
 [ 89  90]
 [ 89 129]]
R at a corner: 1.35e+10   on an edge: -3.51e+09   in a flat area: 0

The detector finds the 8 corners of the two shapes and nothing on the straight edges. Look at the three R values: large positive at a corner, negative on the top edge, and zero in the flat area, exactly as the table predicts.

A close cousin is the Shi-Tomasi detector, which uses the smaller eigenvalue directly as the score. OpenCV's cv2.goodFeaturesToTrack uses it by default and is popular for tracking points in video.

Descriptors: a fingerprint for each keypoint

Finding corners is half the job. To match a corner in photo A with the same corner in photo B, each keypoint needs a descriptor: a short list of numbers that describes the patch of pixels around it. Similar patches should get similar descriptors.

A fast and clever kind is the binary descriptor. Pick, say, 128 random pairs of pixel positions inside the patch. For each pair, write 1 if the first pixel is brighter than the second, else 0. The patch becomes a string of 128 bits. This idea is called BRIEF.

To compare two binary descriptors we count how many bits differ. This count is the Hamming distance: 0 means identical, and around half the bits for unrelated patches. Computers can count differing bits extremely fast.

import numpy as np

rng = np.random.default_rng(0)
PATCH = 15
pairs = rng.integers(0, PATCH, size=(128, 4))          # 128 random (r1, c1, r2, c2) pairs

def describe(img, r, c):
    h = PATCH // 2
    patch = img[r - h:r + h + 1, c - h:c + h + 1]
    return patch[pairs[:, 0], pairs[:, 1]] > patch[pairs[:, 2], pairs[:, 3]]   # 128 bits

def hamming(a, b):
    return int(np.count_nonzero(a != b))

# a textured image, and the same image shifted 5 rows down and 8 columns right, with noise
texture = rng.integers(0, 256, (80, 80)).astype(float)
shifted = np.roll(texture, (5, 8), axis=(0, 1)) + rng.normal(0, 10, texture.shape)

d_original = describe(texture, 30, 30)
d_same_point = describe(shifted, 35, 38)               # the same point after the shift
d_other_point = describe(shifted, 50, 20)              # a different point
print("same point, Hamming distance :", hamming(d_original, d_same_point), "of 128")
print("other point, Hamming distance:", hamming(d_original, d_other_point), "of 128")

# search the whole shifted image for the best match
best = min(((hamming(d_original, describe(shifted, r, c)), r, c)
            for r in range(10, 70) for c in range(10, 70)))
print("best match at row, col:", best[1], best[2], " distance:", best[0])

same point, Hamming distance : 6 of 128
other point, Hamming distance: 59 of 128
best match at row, col: 35 38  distance: 6

The true match has a small Hamming distance despite the noise, while a random point differs in about half the bits. A brute-force search over every position finds the correct shifted location, (35, 38).

ORB: corners plus descriptors, ready to use

ORB (Oriented FAST and Rotated BRIEF) is OpenCV's free, fast feature pipeline:

1. FAST finds corner candidates by checking a circle of 16 pixels around each point: if enough of them in a row are all much brighter or all much darker than the centre, it is a corner. It is very fast.
2. ORB keeps the best candidates (ranked with the Harris score) at several image sizes, so it copes with zoom.
3. It measures each keypoint's orientation (which way the patch "points") and rotates the BRIEF pairs to match. So the descriptor works even if the second photo is rotated.

Here ORB matches a test scene with a copy rotated by 20 degrees. The ratio test keeps a match only if the best candidate is clearly better than the second-best, which removes most wrong matches.

import cv2
import numpy as np

# a textured test scene: random rectangles and circles
rng = np.random.default_rng(1)
scene = np.full((300, 400), 90, np.uint8)
for _ in range(25):
    x, y = rng.integers(0, 360), rng.integers(0, 260)
    cv2.rectangle(scene, (int(x), int(y)), (int(x + rng.integers(15, 40)), int(y + rng.integers(15, 40))),
                  int(rng.integers(0, 256)), -1)
    cv2.circle(scene, (int(rng.integers(0, 400)), int(rng.integers(0, 300))), int(rng.integers(5, 15)),
               int(rng.integers(0, 256)), -1)

# the same scene rotated by 20 degrees (like a second photo from a tilted phone)
M = cv2.getRotationMatrix2D((200, 150), 20, 1.0)
rotated = cv2.warpAffine(scene, M, (400, 300), borderValue=90)

orb = cv2.ORB_create(nfeatures=500)
kp1, des1 = orb.detectAndCompute(scene, None)
kp2, des2 = orb.detectAndCompute(rotated, None)
print("keypoints:", len(kp1), len(kp2), " descriptor shape:", des1.shape, des1.dtype)

bf = cv2.BFMatcher(cv2.NORM_HAMMING)
pairs = bf.knnMatch(des1, des2, k=2)
good = [m for m, n in pairs if m.distance < 0.75 * n.distance]     # ratio test
print("matches after ratio test:", len(good))

# how many good matches agree with the true rotation?
src = np.float32([kp1[m.queryIdx].pt for m in good])
dst = np.float32([kp2[m.trainIdx].pt for m in good])
predicted = cv2.transform(src[None], M)[0]
correct = np.linalg.norm(predicted - dst, axis=1) < 3
print(f"correct matches: {correct.sum()} of {len(good)}")

keypoints: 459 476  descriptor shape: (459, 32) uint8
matches after ratio test: 205
correct matches: 187 of 205

Each ORB descriptor is 32 bytes, which is 256 bits. About 91% of the matches that survive the ratio test are correct. The few wrong ones are normal; in the perspective transform lesson you will see how RANSAC ignores them while estimating how one image maps onto the other.

Where features are used

• Panoramas: match features between overlapping photos, then warp and blend.
• Document and object alignment: line up a scanned form with a template to read fields from fixed positions.
• Augmented reality: find a poster or book cover in the camera view and draw on top of it.
• Visual odometry and SLAM: robots and drones track features frame by frame to estimate how they moved.

Common mistakes

• Matching without a ratio test or RANSAC. Many raw matches are wrong; filter them.
• Using NORM_L2 for ORB. Binary descriptors need cv2.NORM_HAMMING.
• Feature matching on plain, textureless surfaces. A white wall has no corners; features need texture.
• Expecting features to survive huge changes. Big changes in viewpoint, lighting or blur reduce matches sharply.

Key points

• Corners can be located exactly because the image changes in every direction around them.
• Harris uses the structure tensor of gradients; its score R is high at corners, negative on edges.
• A descriptor summarises the patch around a keypoint; binary descriptors compare pixel pairs.
• Hamming distance counts differing bits and is the right distance for ORB and BRIEF.
• ORB = FAST corners + orientation + rotated BRIEF; use cv2.ORB_create, BFMatcher(cv2.NORM_HAMMING) and the ratio test.

import numpy as np
from numpy.lib.stride_tricks import sliding_window_view

def conv2d(img, kernel):
    p = np.pad(img.astype(np.float64), kernel.shape[0] // 2, mode="edge")
    return np.einsum("ijkl,kl->ij", sliding_window_view(p, kernel.shape), kernel)

def gaussian_kernel(size, sigma):
    h = size // 2
    y, x = np.mgrid[-h:h + 1, -h:h + 1]
    k = np.exp(-(x ** 2 + y ** 2) / (2 * sigma ** 2))
    return k / k.sum()

img = np.full((100, 140), 40.0)
img[20:60, 20:70] = 200
img[50:90, 90:130] = 140
rr, cc = np.ogrid[:100, :140]
img[(rr - 25) ** 2 + (cc - 110) ** 2 <= 15 ** 2] = 255

sx = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
gx, gy = conv2d(img, sx), conv2d(img, sx.T)
g = gaussian_kernel(5, 1.0)
Sxx, Syy, Sxy = conv2d(gx * gx, g), conv2d(gy * gy, g), conv2d(gx * gy, g)
R = (Sxx * Syy - Sxy ** 2) - 0.05 * (Sxx + Syy) ** 2
p = np.pad(R, 3, constant_values=-np.inf)
local_max = R == sliding_window_view(p, (7, 7)).max(axis=(2, 3))
print("corners:", len(np.argwhere(local_max & (R > 0.01 * R.max()))))