Harris Corner Detector

The Harris corner detector finds points where image intensity changes strongly in more than one direction. For every pixel it builds a $2 \times 2$ matrix from the image gradients in a surrounding window and computes a corner score from the matrix’s determinant and trace; corners are the local maxima of this score above a threshold. Proposed by Chris Harris and Mike Stephens in 1988, and also known as the Plessey detector after their employer, it is one of the most widely used feature detectors, and its matrix underlies the Shi–Tomasi detector, the Lucas–Kanade method, and the scale- and affine-adapted detectors that followed.

Problem

Find image points that can be localized precisely in both directions, so that they can be found again in another image. In flat regions a small window looks the same wherever it is placed, and along a straight edge it can slide parallel to the edge without changing. Only where the intensity varies in two independent directions, at corners, junctions, and isolated spots, does every shift change the window’s contents.

Inputs and Outputs

• Input: a grayscale image.
• Output: a response map $R(x, y)$ and a list of corner positions, the local maxima of $R$ above a threshold, each with its response value.

Intuition

Place a small window on the image and shift it slightly. Compare the shifted contents with the original by the sum of squared differences. In a flat region the difference stays small for every shift; on an edge it is large for shifts across the edge and small for shifts along it; at a corner it is large for every shift. Moravec’s detector, which Harris and Stephens set out to improve, tested only shifts in four discrete directions with a binary rectangular window and took the minimum. Harris and Stephens replaced the discrete shifts with a Taylor expansion that covers every direction at once, replaced the window with a Gaussian, and replaced the minimum with a score that uses the full shape of the change.

Mathematical Formulation

The structure tensor

The change produced by shifting a window $w$ by $\Delta = (u, v)$ is

$$E(u, v) = \sum_{x, y} w(x, y) \left[ I(x + u, y + v) - I(x, y) \right]^2 .$$

A first-order Taylor expansion, $I(x + u, y + v) \approx I(x, y) + u I_x + v I_y$, turns $E$ into a quadratic form:

$$E(u, v) \approx \begin{bmatrix} u & v \end{bmatrix} M \begin{bmatrix} u \\ v \end{bmatrix}, \qquad 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}.$$

$M$ is the structure tensor, also called the second-moment or auto-correlation matrix. It is the same matrix that appears in the normal equations of window-based motion estimation, and its eigenvalues $\lambda_1 \ge \lambda_2 \ge 0$ classify a window as flat, edge, or corner, as explained for the aperture problem and feature tracking. Harris and Stephens noted that the eigenvalues are proportional to the principal curvatures of the local auto-correlation function.

In modern implementations the derivatives are computed with derivative-of-Gaussian filters at a differentiation scale $\sigma_D$, and the window $w$ is a Gaussian of integration scale $\sigma_I$:

$$M = g_{\sigma_I} * \begin{bmatrix} I_x^2 & I_x I_y \\ I_x I_y & I_y^2 \end{bmatrix}.$$

The 1988 paper used simple $(-1, 0, 1)$ differences for the derivatives.

The corner response

Rather than computing the eigenvalues, Harris and Stephens used the determinant and trace, $\det M = \lambda_1 \lambda_2$ and $\operatorname{trace} M = \lambda_1 + \lambda_2$:

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

$R$ is large and positive when both eigenvalues are large (a corner), negative when one dominates (an edge), and small in magnitude when both are small (a flat region). Writing $r = \lambda_2 / \lambda_1$ gives $R = \lambda_1^2 \left[ r - k (1 + r)^2 \right]$, so $R > 0$ requires the smaller eigenvalue to be at least a certain fraction of the larger: about 4.4% for $k = 0.04$, 5.6% for $k = 0.05$, and 6.9% for $k = 0.06$. Larger $k$ therefore rejects more edge-like points. The paper itself gives no value for $k$; Tuytelaars and Mikolajczyk describe 0.04 as typical, and scikit-image’s corner_harris defaults to 0.05.

Algorithm

1. Compute the derivatives $I_x$ and $I_y$, usually with derivative-of-Gaussian filters at scale $\sigma_D$.
2. Form the three product images $I_x^2$, $I_y^2$, and $I_x I_y$.
3. Smooth each product with a Gaussian of scale $\sigma_I$ to obtain the entries $A$, $B$, $C$ of $M$ at every pixel.
4. Compute $R = (AB - C^2) - k (A + B)^2$.
5. Threshold $R$, either at an absolute value or at a fraction of its maximum.
6. Apply non-maximum suppression: keep pixels that are the maximum of $R$ in their neighborhood. The original paper used the 8-neighborhood.
7. Optionally refine positions to subpixel accuracy and keep the $N$ strongest corners.

Parameters

• $k$: trades corner selectivity against edge rejection, as derived above.
• $\sigma_D$: the scale of the structures whose gradients are measured. Larger values suppress noise but round off fine corners.
• $\sigma_I$: the size of the neighborhood. It should be larger than $\sigma_D$; Mikolajczyk and Schmid set $\sigma_D = 0.7\,\sigma_I$ in Harris–Laplace.
• Threshold: controls the number of corners. Because $R$ scales with the fourth power of contrast (below), a fraction of the maximum response or a fixed number of strongest corners is more portable across images than an absolute value.
• Suppression radius: the minimum spacing between reported corners.

Complexity

All steps are local. With separable Gaussian filters, the cost per pixel is proportional to the filter widths, so a full response map costs $O(N \sigma)$ for an $N$-pixel image, and non-maximum suppression adds $O(N)$. Detectors built for speed, such as FAST, replace the gradient computation with intensity comparisons (see feature detection).

Implementation

A NumPy/SciPy implementation, applied to a noisy synthetic image containing a square and a disk, and to a copy rotated by 30°:

import numpy as np
from scipy import ndimage

def harris_response(img, sigma_d=1.0, sigma_i=2.0, k=0.05):
    # Gaussian derivatives at the differentiation scale sigma_d.
    Ix = ndimage.gaussian_filter(img, sigma_d, order=(0, 1))
    Iy = ndimage.gaussian_filter(img, sigma_d, order=(1, 0))
    # Structure tensor entries: gradient products averaged by a Gaussian window.
    A = ndimage.gaussian_filter(Ix * Ix, sigma_i)
    B = ndimage.gaussian_filter(Iy * Iy, sigma_i)
    C = ndimage.gaussian_filter(Ix * Iy, sigma_i)
    return (A * B - C**2) - k * (A + B)**2

def detect(R, rel_threshold=0.1, nms_size=9):
    # Local maxima of R in a 9x9 window that exceed 10% of the strongest response.
    peaks = (R == ndimage.maximum_filter(R, size=nms_size)) & (R > rel_threshold * R.max())
    return np.argwhere(peaks)  # (row, col)

def label(p, corners, disk_center, disk_radius):
    if np.min(np.linalg.norm(corners - p, axis=1)) < 3:
        return "square corner"
    if abs(np.linalg.norm(p - disk_center) - disk_radius) < 3:
        return "disk edge"
    return "other"

# Synthetic image: a square and a disk on a dark background, plus noise.
rng = np.random.default_rng(0)
n = 160
yy, xx = np.mgrid[0:n, 0:n]
img = np.zeros((n, n))
img[30:70, 30:70] = 1.0
disk_center, disk_radius = np.array([110.0, 105.0]), 22.0
dist_to_center = np.hypot(yy - disk_center[0], xx - disk_center[1])
img[dist_to_center <= disk_radius] = 1.0
img += 0.05 * rng.normal(size=img.shape)
square_corners = np.array([[29.5, 29.5], [29.5, 69.5], [69.5, 29.5], [69.5, 69.5]])

R = harris_response(img)
print("Original image:")
for p in detect(R):
    print(f"  detection at (row={p[0]:3d}, col={p[1]:3d})  R={R[tuple(p)]:.2e}  "
          f"-> {label(p, square_corners, disk_center, disk_radius)}")
ring = np.abs(dist_to_center - disk_radius) < 3
print(f"  max R anywhere near the disk edge: {R[ring].max():.2e}")
for name, p in [("square edge", (50, 30)), ("disk edge", (110, 127)), ("flat area", (15, 140))]:
    print(f"  R at {name:<11} {p}: {R[p]: .2e}")

# Rotate by 30 degrees about the image center and detect again.
angle = 30.0
rot = ndimage.rotate(img, angle, reshape=False, order=1)
t, c = np.deg2rad(angle), (n - 1) / 2
rot_matrix = np.array([[np.cos(t), -np.sin(t)], [np.sin(t), np.cos(t)]])  # (row, col)
rot_corners = (square_corners - c) @ rot_matrix.T + c
rot_disk = (disk_center - c) @ rot_matrix.T + c
print(f"Rotated by {angle:.0f} degrees:")
for p in detect(harris_response(rot)):
    err = np.min(np.linalg.norm(rot_corners - p, axis=1))
    print(f"  detection at (row={p[0]:3d}, col={p[1]:3d})  "
          f"-> {label(p, rot_corners, rot_disk, disk_radius)}, {err:.1f} px from the true corner")

Output:

Original image:
  detection at (row= 31, col= 31)  R=5.53e-04  -> square corner
  detection at (row= 31, col= 68)  R=6.42e-04  -> square corner
  detection at (row= 68, col= 31)  R=6.47e-04  -> square corner
  detection at (row= 68, col= 68)  R=6.97e-04  -> square corner
  max R anywhere near the disk edge: -1.72e-06
  R at square edge (50, 30): -1.45e-04
  R at disk edge   (110, 127): -4.18e-05
  R at flat area   (15, 140):  1.29e-08
Rotated by 30 degrees:
  detection at (row= 43, col= 45)  -> square corner, 2.0 px from the true corner
  detection at (row= 62, col= 13)  -> square corner, 2.0 px from the true corner
  detection at (row= 75, col= 64)  -> square corner, 2.0 px from the true corner
  detection at (row= 94, col= 32)  -> square corner, 2.0 px from the true corner

Exactly four points are detected, one at each corner of the square. The response is negative along the straight edges and everywhere along the disk’s boundary: a gently curved edge still has one dominant gradient direction within the window. The flat background gives a response near zero. After rotation, the same four corners are found at their rotated positions. In both images the maxima lie about two pixels inside the true corner, a localization bias caused by the smoothing at scales $\sigma_D$ and $\sigma_I$.

Properties and Behavior

• Rotation invariance. Rotating the image rotates the gradients and turns $M$ into $Q M Q^\top$ for a rotation matrix $Q$, which has the same eigenvalues, so $R$ is unchanged when the window and derivative filters are isotropic. This is what the analytic expansion and the Gaussian window bought over Moravec’s detector. Discrete filters and resampling make the invariance approximate in practice.
• Partial invariance to intensity changes. Adding a constant to the image leaves the derivatives unchanged. Multiplying the contrast by $a$ multiplies $M$ by $a^2$ and $R$ by $a^4$, so the locations of the maxima are unchanged, but a fixed threshold selects different points. Harris and Stephens noted this contrast scaling when they designed the response.
• No scale invariance. $M$ is computed at fixed scales, so the same structure gives different responses at different zoom levels: a sharp corner viewed from far away and a rounded corner viewed close up look different to the detector. Mikolajczyk and Schmid reported that the repeatability of the standard Harris detector was acceptable only for scale changes up to a factor of 1.4.
• Good localization. Corners are localized in both directions, which is why Harris points are widely used for tracking and calibration, though with the inward bias seen above.

Limitations

• Not affine invariant. Under strong viewpoint change, corners change shape and the circular window covers different scene areas in the two images.
• Threshold dependence. The number of detections depends on contrast, and textured regions produce clusters of corners while plain regions produce none, unless the detections are spread out explicitly.
• No descriptor. The detector gives positions only; matching requires a separate descriptor.

Variants

• Shi–Tomasi. Shi and Tomasi used $\min(\lambda_1, \lambda_2)$ as the score, the condition under which the tracking equations are well conditioned; it removes the parameter $k$ at the cost of computing the eigenvalues. See feature tracking.
• Harmonic mean. The score $2 \det M / \operatorname{trace} M = 2 \lambda_1 \lambda_2 / (\lambda_1 + \lambda_2)$, the harmonic mean of the eigenvalues, with a small constant added to the denominator, needs no $k$. scikit-image offers it as an alternative to $k$ in corner_harris, and its documentation names this form Noble’s corner measure. Tuytelaars and Mikolajczyk note that these alternative corner strengths have been shown to correspond to different matrix norms.
• Multi-scale Harris and Harris–Laplace. Mikolajczyk and Schmid computed a scale-adapted matrix, with the derivatives normalized by $\sigma_D^2$, at a series of scales, detected Harris maxima at each, and kept the points where the scale-normalized Laplacian of Gaussian reached a maximum over scale. The selected scale defines a scale-invariant neighborhood for each point.
• Harris-Affine. The same authors extended Harris–Laplace by iteratively estimating an elliptical neighborhood from the second-moment matrix, giving regions that adapt to affine deformations.
• Edges. The 1988 paper used negative responses to detect edges as well as corners; modern use keeps only the corner part.

For blob-like keypoints with built-in scale selection, see SIFT.