Feature Discovery, Part 3: Harris Corner Discovery

0 0 5 minutes read

Feature Discovery, Part 3: Harris Corner Discovery

Character discovery computer vision domain that focuses on using tools to find regions of interest in images. An important feature of most feature detection algorithms is that they do not use machine learning under the hood, which makes the results even clearer and faster in some cases.

In the previous two articles of this series, we looked at the most popular operators for finding image edges: Sobel, Scharr, The Laplacianas well as Gaussian used for image smoothing. In a way, these operators used under-the-hood derivatives and gradients, which are represented by convolutional kernels.

Like edges, in image analysis, another type of geographic region is often examined: corners. Corners occur more often than edges and usually indicate a change in the direction of an object's boundary or the end of one object and the beginning of another. Corners are rarely found, and they provide very valuable information.

For example

Imagine putting together a 2D puzzle. What most people do initially is to find a slice that contains the part of the image that contains the border (edge) of the object. Why? Because in this way, it is easy to identify adjacent pieces, since the number of pieces that share the same edge is small.

We can go further and focus on selecting not the edges but the corners – the region where the object changes its orientation. These pieces are rarer than rims and allow for an even easier search for other nearby pieces due to their unique shape.

For example, in the puzzle below, there are 6 edges (B2, B3, B4, D2, D3again D4) and only 1 corner (C5). By choosing a corner from the beginning, it is easier to make its place because it is rarer than the edges.

The goal of this article is to understand how corners can be found. To do that, we will understand the details of the algorithm for finding the Harris corner – one of the simplest and most popular methods developed in 1988.

An idea

Let's take three types of states: flat, the edgeagain it exists. We have already shown the structure of these regions above. Our aim will be to understand the distribution of gradients in all three cases.

During our analysis, we will also create an ellipse that contains most of the plotted points. As we will see, its condition will give strong indications of the type of region we are dealing with.

A flat region

A flat surface is an easy case. In general, the entire region of the image has almost the same intensity values, which makes the gradient values along the X and Y axes small and centered on 0.

Taking the gradient points (Gₓ, Gᵧ) from the flat image example above, we can plot their distribution, which looks like this below:

Now we can construct an ellipse at the set points with center at (0, 0). Then we can identify its two main axes:

I the major axis where the ellipse is greatly expanded.
I minor axis where the ellipse reaches its minimum.

In the case of a flat region, it may be difficult to distinguish visually between the major and minor axes, since the ellipse is often circular in shape, as in our case.

However, for each major axis, we can then calculate the ellipse radii λ₁ and λ₂. As shown in the picture above, they are almost equal and have relatively small values.

The edge region

In the edge region, the tension changes only at the edge region. Apart from the edge, the intensity remains almost the same. Given that, most gradient points still lie at (0, 0).

However, in a small area around the threshold, the gradient values can change significantly in both directions. From the example image above, the edge is diagonal, and we can see changes in both directions. Therefore, the gradient distribution is skewed in the diagonal direction as shown below:

In the peripheral regions, the planar ellipse tends to be tilted in one direction and has very different radii λ₁ and λ₂.

A corner district

For corners, most of the stiffness values outside the corners remain the same; therefore, the distribution of most points is still found near the point (0, 0).

If we look at the structure of a corner, we can think of it as the intersection of two edges with two opposite directions. Finally, we have already discussed in the previous section that the distribution goes in the same direction either in X or Y, or both directions.

By having two edges at the corner, we end up with two different spectrums growing in two different directions from the center. An example is shown below.

Finally, if we construct an ellipse around that distribution, we will notice that it is larger than in flat areas or edges. We can isolate this effect by measuring λ₁ and λ₂, which will take much larger values in this case.

To visualize

We have just seen three cases where λ took different values. To better see the results, we can create the diagram below:

Diagram showing the relationship between λ values and circuit types.

The formula

In order to be able to divide a region into one of three zones, the formula below is often used to estimate the R coefficient:

R = λ₁ ⋅ λ₂ – k ⋅ (λ₁ + λ₂)² , where 0.04 ≤ k ≤ 0.06

Based on the R value, we can classify the image region:

R < 0 - edge area
R ~ 0 – flat circuit
R > 0 – corner region

OpenCV

Harris Corner detection can be easily done in OpenCV using the cv2.CornerHarris function. Let's see in the example below how it can be done.

Here's the installation image we'll be working with:

First, let's import the required libraries.

import numpy as np
import cv2
import matplotlib.pyplot as plt

We will convert the input image to grayscale format, as the Harris detector works with pixel density. It is also necessary to convert the image format to float32, since the computational values associated with pixels can exceed the limits. [0, 255].

path = 'data/input/shapes.png'
image = cv2.imread(path)
grayscale_image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
grayscale_image = np.float32(grayscale_image)

Now we can use the Harris filter. I cv2.is thereHarris function takes four parameters:

gray image – input grayscale image in float32 format.
blockSize (= 2) – defines the size of the pixel block in the target pixel area considered for corner detection.
size (= 3) – the magnitude of the Sobel filter used to calculate the derivative.
k (= 0.04) – the coefficient in the formula used to calculate the value of R.

R = cv2.cornerHarris(grayscale_image, 2, 3, 0.04)
R = cv2.dilate(R, None)

The cv2.cornerHarris function returns a matrix of exact dimensions as the original grayscale image. Its values may be outside the normal range [0, 255]. For every pixel, that matrix contains the R value we looked at above.

cv2.dilate is a morphological operator that can quickly select after better grouping local corners.

A common approach is to define a lower limit at which pixels are considered to be corners. For example, we can consider all image pixels as corners whose R value is greater than the global maximum R value divided by 100. In our example, we assign such pixels to the color red (0, 0, 255).

To visualize the image, we need to convert it to RGB format.

image[R > 0.01 * R.max()] = [0, 0, 255]
image_rgb = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)

Finally, we use maplotlib to display the output image.

plt.figure(figsize=(10, 8))
plt.imshow(image_rgb)
plt.title('Harris Corner Detection')
plt.axis('off')
plt.tight_layout()
plt.show()

Here is the result:

Output image. The red color shows the corners.

The conclusion

In this article, we explored a robust method for determining whether an image region exists. The presented formula for estimating the R coefficient works well in most cases.

In real life, there is a common need to implement an edge separator for every image. Building an ellipse around the gradient points and measuring the R coefficient each time is very resource intensive, so more advanced optimization methods are used to speed up the process. However, they are very much based on the theory we read here.