Harris Corner and Edge Detector (2024)

The Harris Corner Detector works well for perfect corners, but what if the corners are not so perfect? Then it might not work so well, but how do we know this? The answer lies in understanding the theory behind the method and recognizing the tradeoffs. In this post we break down the approach and derivation used in the original paper published in 1988, to follow up the derivation we implement the method from scratch to gain a deeper understanding. Code for this tutorial is located on github.

Here’s an overview of the post:

• Background — What is a corner?
• Morvec’s Corner Detector — Where it all began
• Harris Corner Detector — How Harris improved upon Morvec’s
• Implement from Scratch — If you’re here for the code, this one’s for you. Or for the tradeoff.
• Conclusion — tl;dr;

Note: This post assumes that you are familiar with image gradients.

The original intent for the Harris Corner detector was to detect both corners and edges for motion analysis throughout a sequence of images. The reason for tracking edges is that they enable high level features to be obtained since they provide insight to different objects, surfaces, and textures. However, matching edges from frame to frame is a very difficult task, so instead we aim to track corners and fill in with edges.

What exactly is a corner?

In images, a corner is an area where two edges meet together. In fact, the original paper argues that corners should be the corner stone (pun intended) of feature detection because they have the following characteristics:

• Discrete — they take up a unique and small space (not a continuum)
• Reliable — likely to re-appear in a similar area from frame to frame
• Meaningful — strong corners usually come from actual objects rather than shadows or lighting

A corner is an area where two edges meet together

Intuition with Gradients

As a preview to build our intuition, let’s see what the x (horizontal) and y (vertical) gradients of a corner look like. We can compute these with the Sobel operator.

On the leftmost side of figure 1 we have an image of a corner with added noise. The middle two images are the x and y gradients of theimage. The right most side if figure 1 is a scatter plot of the x VS y gradients. The large cluster around the center comes from the flat regions, which have small gradients and are only non-zero due to the added noise. The large values in the x and y gradients create clusters at the right and top respectively, these come from the changes in each direction due to the corner. Let’s take a look at a horizontal edge to better understanding of what’s going on.

In the right most plot of figure 2, we can still see the cluster at the center due to the flat regions, however now we have one smaller cluster above due to the large values in the vertical gradient. The horizontal edge produces a large vertical change (at the edge) as we transverse the image vertically resulting in this cluster on our scatter plot, it is a similar case for the vertical edge. This means that variation between the x and y image gradients can allow us to detect corners.

• High variation in both the x and y gradient directions indicates a corner
• High variation in either the x and y gradient directions indicates an edge
• Low variation in both the x and y gradient directions indicates a flat surface

This is intuitive, but how can we efficiently compute this? The Harris Corner Detector provides a convenient way to detect corners with similar logic.

The Harris Corner Detector is an extension of Morvec’s Corner Detector, which we will now introduce, first intuitively and then Mathematically.

Intuitive Corner Detection with Windows

Consider the perspective of a single window of an image, at the very center of this window we have a pixel under test, we want to determine if this pixel is a corner or not. To approach this, we slide the window in multiple directions and observe the total intensity (sum of all the window pixels) in each direction. In figure 3 below, we can see that there are three different types of image features to consider.

The methods of region classification should feel similar to the gradient intuition above, since image gradients describe the directional change intensity. Morvec’s Method places a window around each pixel and classifies it based on the different changes in intensity as the window is shifted along each direction. However, only four directional shifts are taken into account:

A corner can be declared when the minimum change produced by any of the shifts exceeds a threshold

Mathematical Expression

We can express Morvec’s Method mathematically where a window is centered on (u, v) and is shifted by (x, y). The sum of squared differences (SSD) between all the shifted patch intensities is the total change E.

Drawbacks

1. Only a small set of shifts is considered
2. The response is noisy due the rectangular window
3. Doesn’t consider the entire variation of E, only considers the minimum

The Harris corner detector addresses each of these drawbacks and enables a more robust detection of corners and edges.

Derive an Analytical Expression to consider all possible shifts

To consider every small shift rather than just four shifts, we use a Taylor Expansion about the shift origin (u,v) to derive a convenient analytical expression. This assumes that the horizontal and vertical derivatives of the image exist at (x, y) and (x+u, y+v).

Use a Gaussian Window to reduce noise

The squared first order derivatives are convolved with some window w, in this case we can use a Gaussian window to reduce noise. They are expressed as:

If we use a standard 3x3 window, then there will be 8 total shifts.

Exploit the Quadratic form to consider the total variation

The resulting analytical expression above for the total intensity change E is in a quadratic form and it conveniently allows us to study it’s shape at the origin via eigenvalues of the M matrix.

The eigenvectors are basis vectors that capture the underlying patterns in our data, while the eigenvalues describe the variance in each of the principle (eigenvector) directions. In other words, the eigenvectors scaled by their respective eigenvalues form the axes of an ellipsoid that encompasses a scatter plot of our data. In this case, the eigenvalues describe the change in image intensity in each principle direction. In figure 4 below, we can see that the ellipsoid formed by the eigenvectors/values encompasses the scatter plot of the x VS y image gradients.

The eigenvalues of the M matrix describe the change in image intensity in each of the principle directions

This allows us to assess windowed pixels of an image based on the eigenvalues of the M matrix and we don’t even need the eigenvectors! The magnitudes of the eigenvalues alone are sufficient, since they describe the total variation of gradients (with all directions accounted for)!

The Corner Response contours in figure 5 above show how we can classify each pixel as a flat, edge, or corner (these should look familiar!).

1. Flat — when both eigenvalues are small (E will be smooth)
2. Edge — when one eigenvalue is large and the other is small (E will be ridge shaped)
3. Corner — when both eigenvalues are large (E will be sharply peaked)

Furthermore, we can empirically derive an expression for the corner response that does not involve direct computation of the eigenvalues nor does it involve explicitly considering every (x,y) shift! It only involves computing the Determinant and Trace of matrix M and setting a hyperparameter k that allows us to control which regions are considered flat or not flat.

The Harris corner response R will be positive in corner regions, negative in edge regions, and small in flat regions. Just by using this value, we can find both corners and edges.

The Algorithm to compute the corner response is:

1. Compute the horizontal and vertical image derivatives
2. Compute the product of the derivatives (dx*dx, dy*dy, dx*dy)
3. Convolve a Gaussian Smoothing kernel with each of the product derivatives → Any type of blurring should help reduce noise
4. Compute Corner Response: R = (AB — C²) — k (A + B)²
5. Threshold Corner Response to get strong corners: R > thresh

Let’s Compute the Harris Corner Response on some test images to get a better understanding of what’s going on. The test images are simple edges and corners with a touch of added noise. The code comes from this notebook.

To compute the Corner Response, we will first need the Sobel x and y kernels to compute the image derivatives and a Gaussian Kernel to apply smoothing.

# Sobel kernels
Sx = np.array([
 [1, 0, -1],
 [2, 0, -2],
 [1, 0, -1]])
Sy = Sx.T
# Gaussian Kernel
G = np.array([
 [1, 2, 1],
 [2, 4, 2],
 [1, 2, 1]])/16

Now we can compute the corner response:

def corner_response(image, k, Sx=Sx, Sy=Sy, G=G):
 # compute first derivatives
 dx = cv2.filter2D(image, ddepth=-1, kernel=Sx)
 dy = cv2.filter2D(image, ddepth=-1, kernel=Sy)
 # Gaussian Filter (blur)
 A = cv2.filter2D(dx*dx, ddepth=-1, kernel=G)
 B = cv2.filter2D(dy*dy, ddepth=-1, kernel=G)
 C = cv2.filter2D(dx*dy, ddepth=-1, kernel=G)
 # compute corner response at all pixels
 return (A*B - (C*C)) - k*(A + B)*(A + B)
# hyperparameters
k = 0.05
thresh = 0.5
# thresholded corner responses
strong_corners = corner_response(image, k) > thresh

Figure 7 shows the corner response for each test image without thresholding, notice how the edges are negative while the corners on the bottom have small positive responses at their strongest point. Let’s zoom in and threshold the corners. Notice how more than one pixel is thresholded, we can take care of this by finding the centroid of the response i.e. non-maxima suppression via peak finding (we’ll see how to do this below).

Eigenvalues of the test images

Let’s inspect the eigenvalues of the M matrix at the center pixel of each test image. The image below shows the eigenvectors scaled by their respective eigenvalues for each test image.

Notice that both of the eigenvalues of the flat surface are very small, while the edge eigenvalues each have one large eigenvalue that comes from the change in intensity due to the edge. The eigenvalues of the corner are both relatively large. The blunt corner is very similar to a horizontal edge and may be considered an edge depending on the thresholding, the bluntness really doesn’t make this a good feature so it is not of much worry. The sharp edge has one large eigenvalue and another relatively small eigenvalue that could cause it to be misclassified as an edge, let’s continue to explore this.

Exploring tradeoffs

The sharp edge presents a weakness in the Harris corner detector that is apparent when we look at the eigenvalue and eigenvectors. The rotation of the eigenvectors is relevant to the type of surface (flat, edge, corner), and the corners appear to always have some sort of rotation. This sharp edge seems like a good feature to track, but the Harris Detector may misclassify it due to one of it’s eigenvalues being small! Remember that computing the eigenvectors is computationally expensive, and the Harris Detector cleverly avoids it. This tradeoff balances accuracy with computational efficiency,butmoreimportantly understanding it opens thedoor tostudyingnewmethods.

Full example

Now here’s a full example where we will place points on the corners of a real image. This example is implemented in this notebook.

def get_harris_corners(image, k=k,Sx=Sx, Sy=Sy, G=G):
 # compute corner response
 R = corner_response(image, k, Sx=Sx, Sy=Sy, G=G)
 # find centroids
 ret, labels, stats, centroids = cv2.connectedComponentsWithStats(np.uint8(R > 1e-2))
 # define the criteria to stop and refine the corners
 criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 100, 0.001)
 return cv2.cornerSubPix(image, np.float32(centroids), (9,9), (-1,-1), criteria)

from skimage.data import brick
image = brick().astype(np.float32)
# (OPTIONAL) blurr image --> This will change hyperparameters
# image = cv2.blur(image, ksize=(5,5))
# 0-1 normalize
image /= image.max()
# find corners
corners = get_harris_corners(image)
# draw corners on output image
image_out = np.dstack((image, image, image))
for (x, y) in corners:
 x = np.round(x).astype(int)
 y = np.round(y).astype(int)
 cv2.circle(image_out, (x, y), radius=3, color=(0, 255, 0), thickness=-1)

We can also compute the Harris Edge response by considering a negative threshold:

R = corner_response(image, k)
edge_response = R < -0.025

In this post we learned how the Harris Corner Response can conveniently be computed from products of image derivatives. The corner response can be thresholded to detect corners (positive values) or edges (negative values). We have learned that the derived Corner Response is based on the eigenvalues of the M matrix, where the eigenvalues represent the variance in each gradient direction. We explored a tradeoff, where the method may not be able to detect distinct features such as very sharp edges. Lastly we have seen how to centroid and detect strong corner pixels using opencv in order to display corners as points.

Next Steps

Unpack the corner_response() function and see what the derivatives look like for each test image. Find out what happens when you don’t apply the Gaussian Kernel to the derivative products. Explore some more test images, see what their eigenvalues look like and try to find and understand cases where the method flourishes and cases where it fails.

[1] Harris, C., & Stephens, M. (1988). A combined corner and Edge Detector. Procedings of the Alvey Vision Conference 1988. https://doi.org/10.5244/c.2.23

[2] Peressini, A. L., Sullivan, F. E., & Uhl, J. J. (1993). The mathematics of Nonlinear Programming. Springer.

Thanks for Reading! If you found this useful please consider clapping 👏