Corners, Blobs, and Scale: Finding Repeatable Keypoints

What Makes a Good Keypoint?

In the previous guide we learned to find edges — the places where image brightness changes sharply — using gradients. Edges are a great start, but they are not enough to anchor one photo to another. Imagine you take two pictures of the same building from slightly different spots and want a computer to figure out which point in photo A is the same physical point in photo B. An edge running along a roofline looks identical at every position along that line, so the computer cannot tell where on the edge it is. We need something more pinpointable: a keypoint.

A good image keypoint is judged by four design goals. Repeatability: if the same scene is photographed again — from another angle, in different light, at another zoom — the detector should re-find the very same physical points. Distinctiveness: the little patch around the keypoint should look locally unusual, so it can later be told apart from other patches. Locality: the keypoint should depend only on a small region of the image, so that if part of the scene is hidden (occluded), most keypoints elsewhere survive. Invariance: the keypoint and its measured location should stay stable under rotation, brightness changes, and scale changes. Hold on to these four words — every detector in this guide is an attempt to satisfy them.

Recall the aperture problem from the last guide: peering at an edge through a small window, you can tell how far it moved across the edge, but not how far it slid along it — the view looks unchanged. That is exactly why an edge is poorly localized in 2-D: it pins down only one direction. A corner, where two edges meet, is different: nudge the window in any direction and the patch changes. Because intensity varies in every direction, a corner is locked down in both x and y. This single observation is the seed of the harris corner detector we build next.

The Harris Corner Detector: Intuition

Let's build the corner detector from a single question, no equations yet. Place a small square window on the image. Now slide that window a tiny amount in some direction and ask: *how much did the brightness pattern inside the window change?* If the patch barely changes, the window was sitting on something featureless; if it changes a lot, the window was sitting on something with structure. The trick of Harris is to ask this question for every sliding direction at once and look for places that change a lot no matter which way you go.

Flat region (e.g. a clear patch of sky): sliding the window in any direction barely changes the patch — the response is small everywhere. Not a keypoint.
Edge (e.g. the boundary of a roof against the sky): sliding along the edge barely changes the patch, but sliding across it changes a lot. Large change in one direction, tiny in another. Still ambiguous along the edge.
Corner (e.g. the tip where two roof edges meet): sliding the window in every direction changes the patch substantially. Large change in all directions — this is the keypoint we want.

To turn that thought experiment into a number, we measure the total squared change in the patch when we shift the window by a small amount (u, v). Squaring makes every change count as positive (a pixel that gets darker and one that gets brighter both contribute), and summing over the window aggregates the whole patch into one score. Call this score E(u, v): the energy of the change for shift (u, v).

E(u,v)=\sum_{x,y} w(x,y)\,\bigl[\,I(x+u,\;y+v)-I(x,y)\,\bigr]^2

Let's read this term by term. The sum ∑ runs over every pixel (x, y) inside the window. I(x, y) is the brightness at pixel (x, y) in the original position; I(x+u, y+v) is the brightness at the same patch pixel after the window has shifted by (u, v). Their difference is how much that one pixel's brightness changed under the shift; squaring it makes it a positive 'amount of change'. w(x, y) is a window weight — often a Gaussian bump that counts the center of the window more than the rim, so a single noisy edge pixel can't dominate. So E(u, v) is the weighted total change for one candidate shift. A corner is a location where E(u, v) is large for every direction of (u, v); a flat region has E ≈ 0 everywhere; an edge has E large for shifts across it and ≈ 0 along it. Computing E for all shifts is expensive, so next we use a first-order Taylor approximation — I(x+u, y+v) ≈ I(x, y) + u·I_x + v·I_y — to replace the shifted image with the gradients we already know how to compute. That single substitution turns E into a compact matrix and gives us the harris corner detector.

The Structure Tensor and Corner Response

Plug the Taylor approximation into E. The bracket [I(x+u, y+v) − I(x, y)] becomes simply (u·I_x + v·I_y), the change predicted by the local gradient. Squaring that gives u²I_x² + 2uv·I_x I_y + v²I_y². Every term is a product of gradients times a power of u or v, so when we sum over the window we can factor the (u, v) part out and collect all the gradient sums into a single 2×2 matrix M. The result is wonderfully tidy: E(u, v) ≈ [u v] M [u v]ᵀ. All the information about how the patch changes in every direction now lives inside M, called the second-moment matrix or structure tensor.

M=\sum_{x,y} w(x,y)\begin{bmatrix} I_x^2 & I_x I_y \\ I_x I_y & I_y^2 \end{bmatrix}

Here I_x and I_y are the horizontal and vertical image gradients from the last guide (how fast brightness rises going right, and going down). For each pixel we form three products — I_x² (horizontal edge strength), I_y² (vertical edge strength), and I_x I_y (how the two correlate) — and w(x, y) weights them as we sum over the window. M is symmetric, so it has two real eigenvalues λ₁ and λ₂ with perpendicular eigenvectors. The eigenvectors point along the two principal directions of the patch, and each eigenvalue measures how strongly the patch's brightness varies along its direction. That is the geometric heart of Harris: λ₁ and λ₂ are the 'change in direction 1' and 'change in direction 2'. Two large eigenvalues mean the patch changes a lot in two independent directions — a corner. One large, one small means change in only one direction — an edge. Both small means a flat region.

R=\det(M)-k\,(\operatorname{trace} M)^2=\lambda_1\lambda_2-k(\lambda_1+\lambda_2)^2

Computing eigenvalues at every pixel is slow, so Harris uses a clever shortcut that needs only sums we already have. Two facts from linear algebra: det(M) = λ₁λ₂ and trace(M) = λ₁ + λ₂ — and crucially, det and trace can be read straight off M's entries (det = I_x²·I_y² − (I_x I_y)², trace = I_x² + I_y²) without ever solving for the eigenvalues. The corner response R combines them, with k an empirical constant (typically 0.04–0.06) that tunes how harshly edges are penalized. Read R like a thermometer: large positive R ⇒ corner, large negative R ⇒ edge, |R| small ⇒ flat. Let's classify three patches with k = 0.05. *Flat:* λ₁ = 0.01, λ₂ = 0.02 → R = 0.0002 − 0.05·(0.03)² ≈ +0.00016, tiny → flat. *Edge:* λ₁ = 5, λ₂ = 0.01 → R = 0.05 − 0.05·(5.01)² ≈ −1.2, large negative → edge. *Corner:* λ₁ = 5, λ₂ = 4 → R = 20 − 0.05·(9)² = 20 − 4.05 = 15.95, large positive → corner. Finally we keep only pixels whose R exceeds a threshold and are a local maximum of R (non-maximum suppression), so each corner yields one crisp point instead of a fat blob. That whole pipeline is the harris corner detector.

# Harris corner response (no eigenvalues needed)
import numpy as np
from scipy.ndimage import gaussian_filter, sobel

def harris_response(img, k=0.05, sigma=1.0):
    Ix = sobel(img, axis=1)              # horizontal gradient I_x
    Iy = sobel(img, axis=0)              # vertical gradient I_y
    # Window-weighted sums of gradient products (the Gaussian plays the role of w)
    Sxx = gaussian_filter(Ix * Ix, sigma)
    Syy = gaussian_filter(Iy * Iy, sigma)
    Sxy = gaussian_filter(Ix * Iy, sigma)
    det_M   = Sxx * Syy - Sxy * Sxy      # det(M) = lambda1 * lambda2
    trace_M = Sxx + Syy                  # trace(M) = lambda1 + lambda2
    return det_M - k * trace_M**2        # R: large + at corners, large - at edges

Harris corner response computed from gradient products and a Gaussian window — never solving for eigenvalues directly.

Blobs: Features with Size

Corners are wonderful, but the world is full of features that are not corners. A printed dot, a distant streetlight at night, a cell under a microscope, a coffee stain — these are blobs: roughly circular regions that are uniformly brighter (or darker) than their surroundings. Blob detection looks for exactly these compact bumps of intensity. Unlike a corner, a blob has no sharp meeting of edges; its defining trait is a region that stands out from its background.

Two reasons make blobs essential. First, many natural and man-made structures simply are blob-like, and a corner detector would miss them or fire only weakly on their rims. Second — and this is the idea that powers the rest of the guide — a blob has an intrinsic size. A corner is a zero-dimensional point with no natural width, but a blob is small or large, and that size is a real, measurable property of the feature. Detecting a blob therefore also tells you how big it is, which is the bridge to the idea of scale.

How do we make a filter that lights up on a blob of a given size? Use the Laplacian of Gaussian (LoG). First blur the image a little with a Gaussian (this sets the size the filter cares about), then apply the Laplacian, a second-derivative operator that responds to how much a pixel differs from the average of its neighbours. The combined filter has a 'Mexican-hat' shape: positive in the middle, negative in a surrounding ring. Slide it over the image and it acts like a center-surround detector — it produces a strong response exactly where a bright center is hugged by a darker surround, i.e. on a blob. And it 'rings loudest' when the hat's width matches the blob's width.

\sigma^2\,\nabla^2 G * I,\qquad \nabla^2 G=\frac{\partial^2 G}{\partial x^2}+\frac{\partial^2 G}{\partial y^2}

Reading it: G is a Gaussian blur of width σ (sigma), the knob that sets the filter's size; ∇²G is its Laplacian — the sum of the second derivatives in x and in y — which is the Mexican-hat shape; the \* is convolution, meaning we slide that hat over the whole image I and record its response at every pixel. The leading σ² is the all-important scale normalization. Without it, the Laplacian's response naturally shrinks as you blur more, so big blobs would always look weaker than small ones; multiplying by σ² rescales the responses so they are comparable across sizes. The payoff: as you try the detector at many sizes σ, the response at a true blob's centre peaks at the σ that matches the blob's actual radius. So the LoG doesn't just find blobs — it finds their size. (Concretely, a bright disk of radius r gives the strongest scale-normalized response at about σ ≈ r/√2.)

Scale Space: Seeing at Every Zoom Level

We just saw that a blob's response depends on the filter's size σ. But there is a deeper problem: the same object can fill the frame in one photo and be a tiny speck in another, depending on how far away the camera is. A detector that only works at one size would find a face up close and miss the identical face across the room. The fix is a scale space: instead of analysing the image at a single resolution, we analyse it at a whole continuum of blur levels, and let the data tell us at which scale each feature lives.

To build a scale space, take the image and blur it with Gaussians of increasing width σ, producing a stack of progressively softer copies — a Gaussian pyramid. Increasing σ is like stepping back from the picture: fine detail melts away and only the larger structures survive. A small σ keeps the texture of individual bricks; a large σ keeps only the overall shape of the wall. By stacking many σ values we get a 3-D volume L(x, y, σ): width, height, and now scale as a third axis. A feature that is small in the image shows up at small σ; a big feature shows up at large σ.

\begin{aligned} L(x,y,\sigma) &= G(x,y,\sigma) * I(x,y) \\ D(x,y,\sigma) &= L(x,y,k\sigma) - L(x,y,\sigma) \end{aligned}

The first line defines the scale space: L(x, y, σ) is the image I convolved (\* , slid-and-summed) with a Gaussian G of width σ — i.e. the image blurred by amount σ. The second line is the money-saver. Computing the scale-normalized LoG at many scales is costly, but it turns out the LoG is closely approximated by simply subtracting two neighbouring blur levels: the Difference of Gaussians (DoG), D(x, y, σ) = L(x, y, kσ) − L(x, y, σ), where k is the constant multiplicative step between adjacent scales (e.g. k ≈ 1.26, so each level is about 26 % blurrier than the last). Subtraction is almost free, and the difference of two Gaussians has the same Mexican-hat shape as the LoG, so D is a cheap, accurate blob detector that already carries the σ²-normalization baked in. A strong DoG value means 'a blob of this size lives here.'

Finding Keypoints Across Scales

Now we assemble everything into one detector. We have the DoG volume D(x, y, σ): a stack of difference images, one per scale, sitting in a 3-D space of (x, y, scale). A keypoint is declared at any voxel that is a local extremum — a maximum or minimum — of D in both space and scale at once. Concretely, compare each candidate pixel against its 8 neighbours in the same scale plus the 9 pixels directly above and the 9 directly below in the adjacent scales: 8 + 9 + 9 = 26 neighbours. If the centre value is bigger (or smaller) than all 26, it is a blob that is best-localized both in position and in size. This is how blob detection and scale space combine to pin a feature down in all three coordinates.

# A keypoint is a local extremum of the DoG in space AND scale
def is_keypoint(D, x, y, s, contrast_thresh=0.03):
    # D is the DoG stack indexed [scale, row, col]
    block = D[s-1:s+2, y-1:y+2, x-1:x+2]   # 3 x 3 x 3 = 27 values
    center = D[s, y, x]
    # Compare the center against all 26 neighbours
    if center < block.max() and center > block.min():
        return False                        # neither a max nor a min -> reject
    if abs(center) < contrast_thresh:
        return False                        # too faint -> unstable -> reject
    return True                             # accept candidate (x, y, scale = s)

A DoG keypoint is a value larger or smaller than all 26 of its space-and-scale neighbours, and not too faint.

Raw extrema are good candidates, but master-grade detectors refine them in three plain steps. Sub-pixel / sub-scale localization: fit a small quadratic (a parabola-like surface) to D around the extremum and jump to its true peak, so the keypoint lands between pixels and between scales rather than snapping to the grid. Low-contrast rejection: throw away points whose refined DoG value is too small — they are faint and easily flipped by noise, hurting repeatability. Edge rejection: the DoG also fires along edges, which are poorly localized (remember the aperture problem), so we discard them using the ratio of the two eigenvalues of the local Hessian (the 2×2 matrix of second derivatives). When one eigenvalue hugely outweighs the other, we're on an edge, not a blob — this is the very same eigenvalue-ratio idea Harris used to separate corners from edges. The surviving points form a stable keypoint set.

The output of all this is precisely what the next stage needs: each keypoint is a tuple (x, y, scale) — a location and a characteristic size. The scale matters enormously: it tells the descriptor exactly how big a region around the point to look at, so that a feature measured up close and the same feature measured far away end up described over proportionally sized patches and therefore match. In the next guide we feed each (x, y, scale) keypoint to the SIFT descriptor, which also assigns an orientation and then summarizes the local gradients into a vector you can compare against keypoints from another image.