Harris Corner Detection

Consider a small image patch centered at pixel $(x,y)$, and shift it by a displacement $(u,v)$. The auto-correlation (self-similarity) of the image under this shift is measured by the sum of squared differences between the original and shifted patch:

where $W$ is a local window and $w$ is a weighting function (often a Gaussian, giving more importance to pixels near the center). If $E(u,v)$ is small for all $(u,v)$ in some neighborhood of the origin, the patch is flat or nearly so under shifts — not a corner. If $E$ grows quickly in every direction, the point is a corner. To make this tractable, Taylor-expand the shifted image to first order:

where $I_x, I_y$ are the horizontal and vertical image gradients. Substituting into $E$:

where $M$ is the second-moment matrix (structure tensor):

Here $\langle\cdot\rangle$ denotes the weighted local average over the window. $M$ is a real symmetric $2\times2$ matrix, so it has an eigendecomposition $M = R^{-1}\mathrm{diag}(\lambda_1,\lambda_2)R$ with real eigenvalues $\lambda_1 \ge \lambda_2 \ge 0$ (positive semi-definite because $M$ is a sum of outer products). Geometrically, $E(u,v) \approx \lambda_1 u'^2 + \lambda_2 v'^2$ in the rotated coordinate frame aligned with $M$'s eigenvectors — an elliptical bowl whose axis lengths are governed by $1/\sqrt{\lambda_1}, 1/\sqrt{\lambda_2}$. This yields a clean geometric classification:

• $\lambda_1 \approx \lambda_2 \approx 0$: flat region — the surface $E$ is nearly zero everywhere; no direction changes the patch.
• $\lambda_1 \gg \lambda_2 \approx 0$: edge — the bowl is a trough, flat along the edge direction (eigenvector for $\lambda_2$), steep across it.
• $\lambda_1 \approx \lambda_2 \gg 0$: corner — the bowl is steep in every direction, so any shift increases $E$ substantially.

Computing eigenvalues at every pixel is possible but relatively expensive and requires care with numerical edge cases. Harris and Stephens's key trick is to score "cornerness" using only the trace and determinant of $M$, both of which are cheap, closed-form functions of its entries and invariant under rotation:

The Harris corner response is defined as:

with $k$ an empirical sensitivity constant, typically $k \in [0.04, 0.06]$. Expanding in terms of eigenvalues:

• If both eigenvalues are small, both terms are small, so $R \approx 0$ (flat).
• If one eigenvalue dominates ($\lambda_1 \gg \lambda_2 \approx 0$), then $\lambda_1\lambda_2 \approx 0$ but $(\lambda_1+\lambda_2)^2 \approx \lambda_1^2$ is large, so $R$ becomes strongly negative (edge).
• If both eigenvalues are large and comparable, $\lambda_1\lambda_2$ dominates over $k(\lambda_1+\lambda_2)^2$ (since $k < 1/4$ suppresses the square term relative to the product for balanced eigenvalues), so $R$ is positive and large (corner).

This gives a single scalar per pixel: large positive $R$ signals a corner, large negative $R$ signals an edge, and small $|R|$ signals a flat region — without ever diagonalizing $M$.

Full pipeline: 1. Compute image gradients $I_x, I_y$ (e.g., via Sobel filters). 2. Form the products $I_x^2$, $I_y^2$, $I_xI_y$ at every pixel. 3. Smooth each product with a Gaussian (or box) window to get the windowed sums $\langle I_x^2\rangle, \langle I_y^2\rangle, \langle I_xI_y\rangle$ — this is the local averaging that assembles $M(x,y)$ per pixel. 4. Compute $R(x,y) = \det M - k\,\mathrm{trace}(M)^2$ at every pixel. 5. Threshold $R$ to keep only strong positive responses. 6. Apply non-maximum suppression so each physical corner yields one detected point rather than a cluster. ``` function HARRIS_CORNERS(I, window_sigma, k, threshold, nms_radius):

# Step 1: gradients Ix, Iy <- SOBEL_GRADIENTS(I)

# Step 2: gradient products Ixx <- Ix Ix Iyy <- Iy Iy Ixy <- Ix * Iy

# Step 3: windowed sums (structure tensor entries), per pixel Sxx <- GAUSSIAN_BLUR(Ixx, window_sigma) Syy <- GAUSSIAN_BLUR(Iyy, window_sigma) Sxy <- GAUSSIAN_BLUR(Ixy, window_sigma) R <- array of zeros, same shape as I

# Step 4: response function, per pixel for each pixel (x, y): det_M <- Sxx[x,y] Syy[x,y] - Sxy[x,y]^2 trace_M <- Sxx[x,y] + Syy[x,y] R[x,y] <- det_M - k trace_M^2

# Step 5: threshold candidates <- { (x,y) : R[x,y] > threshold }

# Step 6: non-maximum suppression corners <- [] for (x, y) in candidates, sorted by R[x,y] descending: if R[x,y] is the maximum value of R within the (2nms_radius+1)x(2nms_radius+1) neighborhood of (x,y): corners.append((x, y, R[x,y])) return corners ``` Non-maximum suppression matters because $R$ varies smoothly near a true corner, producing many nearby pixels above threshold; without suppression, the detector would report a blob of points instead of one location per corner.

Time complexity. Let the image have $N = W \times H$ pixels and let the smoothing/gradient kernels have fixed size (independent of $N$, e.g. a $5\times5$ Gaussian). Every stage — gradient computation, forming products, Gaussian smoothing, computing $R$, thresholding — touches each pixel a constant number of times with a constant-size convolution kernel. Hence the entire pipeline is:

If the smoothing kernel is separable (Gaussians are) and has spatial extent $k$, exact convolution costs $O(Nk)$ using two 1D passes rather than $O(Nk^2)$ for a naive 2D convolution — a standard and important constant-factor optimization, but it does not change the asymptotic linearity in $N$. Non-maximum suppression over a window of radius $r$ costs $O(Nr^2)$ naively, or $O(N)$ amortized with a sliding-window maximum algorithm (e.g., using a monotonic deque per row and per column), since $r$ is a fixed small constant in practice.

Space complexity. All intermediate arrays ($I_x, I_y, I_x^2, I_y^2, I_xI_y, S_{xx}, S_{yy}, S_{xy}, R$) are the same size as the image, so space is $O(N)$; in an optimized streaming implementation, several of these can be computed and discarded row-by-row, but the asymptotic bound remains $O(N)$.

Correctness argument. The core claim is that $E(u,v) \approx u^\top M u$ is a faithful local model of auto-correlation, and that the sign/magnitude pattern of $(\lambda_1, \lambda_2)$ correctly separates flat/edge/corner regions. Three facts underwrite this: 1. Positive semi-definiteness of $M$. $M$ is a weighted sum of rank-1 outer products $\begin{pmatrix}I_x\\I_y\end{pmatrix}\begin{pmatrix}I_x&I_y\end{pmatrix}$ with nonnegative weights $w \geq 0$. A nonnegative sum of PSD matrices is PSD, so $\lambda_1,\lambda_2 \geq 0$ always — the model never predicts a decreasing $E$, consistent with $E$ being a sum of squares. 2. Rotation invariance. $\det(M)$ and $\mathrm{trace}(M)$ are invariant under any orthogonal change of coordinates (rotating the image), since $\det$ and trace are similarity-invariants: $\det(P^\top M P) = \det(M)$ and $\mathrm{trace}(P^\top M P) = \mathrm{trace}(M)$ for orthogonal $P$. Therefore $R$, and hence the corner classification, does not depend on the orientation of the image or the local gradient directions — a genuine rotational invariance, not merely an approximate one. 3. $R$ correctly orders the three regimes for $k < 1/4$. Write $R = \lambda_1\lambda_2 - k(\lambda_1+\lambda_2)^2$. Substituting the ratio $\rho = \lambda_2/\lambda_1 \in [0,1]$ (assuming $\lambda_1 \geq \lambda_2 \geq 0$, $\lambda_1 > 0$):

The bracket is positive when $\rho$ is close to 1 (comparable eigenvalues — corner) and negative when $\rho \to 0$ (one eigenvalue vanishes — edge), for any $k \in (0, 1/4)$, because at $\rho=1$ the bracket equals $1 - 4k > 0$, and at $\rho = 0$ it equals $-k < 0$, with the bracket monotonically increasing in $\rho$ on $[0,1]$ (its derivative is $1 - 2k(1+\rho) > 1 - 4k > 0$). This proves the sign of $R$ transitions exactly once as $\rho$ increases from 0 to 1, validating $R$'s use as a corner/edge discriminator for the standard choice of $k$.

Honest tradeoffs. The first-order Taylor expansion of $E(u,v)$ is only locally accurate; for large shifts or highly curved intensity surfaces it degrades, but Harris only needs the local quadratic shape near $u=v=0$, so this is not a practical limitation. The constant $k$ is empirical, not derived from first principles, and different values shift the edge/corner boundary. Gaussian window size trades off localization (small window, precise but noise-sensitive) against robustness (large window, smooth but blurs nearby corners together). Most importantly, the detector is invariant to rotation and (to a good approximation) to affine intensity changes, but not scale-invariant: a corner detected at one image resolution may vanish or shift at another, motivating scale-space treatments such as SIFT Scale-Invariant Features.

Consider a synthetic $5\times 5$ grayscale patch containing a clean right-angle corner (bright region in the top-right, dark elsewhere), centered on the pixel of interest:

The corner of the bright region sits near pixel $(3,2)$ (0-indexed columns, rows). Compute simple forward-difference gradients $I_x = I(x+1,y) - I(x,y)$, $I_y = I(x,y+1)-I(x,y)$ over a $3\times3$ window centered at $(3,2)$, using a uniform (unweighted) window for simplicity:

Sum the columns to get the (unweighted) second-moment matrix entries:

Compute $\det(M)$ and $\mathrm{trace}(M)$ directly (no eigendecomposition needed):

With $k = 0.04$:

$R > 0$ and large — the pixel is correctly classified as a corner. For comparison, repeat the same computation centered at pixel $(3,1)$, which sits on the vertical edge between dark and bright regions but away from the corner tip. There, gradients are strong horizontally ($I_x \approx 100$) but $I_y \approx 0$ throughout the window, giving roughly $M \approx \mathrm{diag}(20000, 0)$, so $\det(M) \approx 0$, $\mathrm{trace}(M) \approx 20000$, and:

Strongly negative $R$ — correctly classified as an edge, exactly as the theory predicts.

The Harris detector emerged from British computer vision research in the late 1980s, building on Hans Moravec's earlier "cornerness" measure, which used four discrete shift directions and suffered from anisotropy and noise sensitivity. Harris and Stephens generalized Moravec's idea to continuous shifts via the Taylor expansion and the second-moment matrix, producing an operator that is smooth, rotation-invariant, and robust — properties Moravec's discrete version lacked.

Harris corners (and its direct descendants, the Shi–Tomasi "minimum eigenvalue" criterion and Harris-Laplace scale-space variants) are foundational building blocks across classical and modern vision pipelines. Camera calibration routines use them to locate checkerboard corners precisely. Structure-from-motion and visual SLAM systems use Harris-style corner scores to select trackable points before running Lucas-Kanade Optical Flow or feature matching across frames, since a good track requires a point that is locally distinctive in both image dimensions. Panorama stitching software historically detected Harris corners as candidate keypoints before computing rotation- and scale-invariant descriptors around them. Optical flow and KLT (Kanade–Lucas–Tomasi) trackers explicitly select points to track using the eigenvalues of the same second-moment matrix, since a well-conditioned $M$ (both eigenvalues large) is precisely the condition under which the optical flow equations are numerically stable — the corner criterion and the "good point to track" criterion are, mathematically, the same object viewed from two purposes. Even in the deep learning era, learned keypoint detectors (e.g., SuperPoint-style networks) are frequently evaluated against Harris corners as a classical baseline, and the underlying structure-tensor idea persists in edge-aware image filters and anisotropic diffusion.

Ask the tutor: How does the Harris corner response change if the image is rotated by 45° — and why does the math guarantee the same corners are detected?