Harris Corner Detector in Python

Thought I'd share a simple Python implementation of the Harris corner detector. I have seen people looking for a python implementation for a range of applications so I'm hoping someone finds this useful.

The Harris (or Harris & Stephens) corner detection algorithm is one of the simplest corner indicators available. The general idea is to locate points where the surrounding neighborhood shows edges in more than one direction, these are then corners or interest points. The algorithm is explained here. In short, a matrix W is created from the outer product of the image gradient, this matrix is averaged over a region and then a corner response function is defined as the ratio of the determinant to the trace of W.

Let's see what this looks like in code. First we need to be able to do convolutions of 2D signals. For this NumPy is not enough and we need to use the signal module in SciPy. Create a file filtertools.py and add the following functions needed to create Gaussian derivative kernels and apply them to the image.

from scipy import *
from scipy import signal


def gauss_derivative_kernels(size, sizey=None):
    """ returns x and y derivatives of a 2D 
        gauss kernel array for convolutions """
    size = int(size)
    if not sizey:
        sizey = size
    else:
        sizey = int(sizey)
    y, x = mgrid[-size:size+1, -sizey:sizey+1]

    #x and y derivatives of a 2D gaussian with standard dev half of size
    # (ignore scale factor)
    gx = - x * exp(-(x**2/float((0.5*size)**2)+y**2/float((0.5*sizey)**2))) 
    gy = - y * exp(-(x**2/float((0.5*size)**2)+y**2/float((0.5*sizey)**2))) 

    return gx,gy

def gauss_derivatives(im, n, ny=None):
    """ returns x and y derivatives of an image using gaussian 
        derivative filters of size n. The optional argument 
        ny allows for a different size in the y direction."""

    gx,gy = gauss_derivative_kernels(n, sizey=ny)

    imx = signal.convolve(im,gx, mode='same')
    imy = signal.convolve(im,gy, mode='same')

    return imx,imy

The point of using Gaussian derivative filters is that this computes a smoothing of the image, to a scale defined by the size of the filter, and the derivatives at the same time. The derivatives are less noisy than if computed with a simple difference filter on the original image.

First add the corner response function to a file harris.py which will make use of the Gaussian derivatives above.

def compute_harris_response(image):
    """ compute the Harris corner detector response function 
        for each pixel in the image"""

    #derivatives
    imx,imy = filtertools.gauss_derivatives(image, 3)

    #kernel for blurring
    gauss = filtertools.gauss_kernel(3)

    #compute components of the structure tensor
    Wxx = signal.convolve(imx*imx,gauss, mode='same')
    Wxy = signal.convolve(imx*imy,gauss, mode='same')
    Wyy = signal.convolve(imy*imy,gauss, mode='same')

    #determinant and trace
    Wdet = Wxx*Wyy - Wxy**2
    Wtr = Wxx + Wyy

    return Wdet / Wtr

This gives an image with each pixel containing the value of the Harris response function. Now it is just a matter of picking out the information needed from this image. Picking all values above a threshold with the additional constraint that corners must be separated with a minimum distance is an approach that often gives good results. To do this, take all candidate pixels, sort them in descending order of corner response values and mark off regions too close to positions already marked as corners. Add this function to harris.py.

def get_harris_points(harrisim, min_distance=10, threshold=0.1):
    """ return corners from a Harris response image
        min_distance is the minimum nbr of pixels separating 
        corners and image boundary"""

    #find top corner candidates above a threshold
    corner_threshold = max(harrisim.ravel()) * threshold
    harrisim_t = (harrisim > corner_threshold) * 1

    #get coordinates of candidates
    candidates = harrisim_t.nonzero()
    coords = [ (candidates[0][c],candidates[1][c]) for c in range(len(candidates[0]))]
    #...and their values
    candidate_values = [harrisim[c[0]][c[1]] for c in coords]

    #sort candidates
    index = argsort(candidate_values)

    #store allowed point locations in array
    allowed_locations = zeros(harrisim.shape)
    allowed_locations[min_distance:-min_distance,min_distance:-min_distance] = 1

    #select the best points taking min_distance into account
    filtered_coords = []
    for i in index:
        if allowed_locations[coords[i][0]][coords[i][1]] == 1:
            filtered_coords.append(coords[i])
            allowed_locations[(coords[i][0]-min_distance):(coords[i][0]+min_distance),(coords[i][1]-min_distance):(coords[i][1]+min_distance)] = 0

    return filtered_coords

Now you have all you need to detect corner points in images. To make it easier to show the corner points in the image you can add a plotting function using matplotlib (PyLab) as follows.

def plot_harris_points(image, filtered_coords):
    """ plots corners found in image"""
    from pylab import *

    figure()
    gray()
    imshow(image)
    plot([p[1] for p in filtered_coords],[p[0] for p in filtered_coords],'*')
    axis('off')
    show()

Try running the following commands on an example image:

>>> im = array(Image.open('empire.jpg').convert("L"))
>>> harrisim = harris.compute_harris_response(im)
>>> filtered_coords = harris.get_harris_points(harrisim,6)
>>> harris.plot_harris_points(im, filtered_coords)

The image is opened and converted to grayscale. Then the response function is computed and points selected based on the response values. Finally, the points are plotted overlaid on the original image. This should give you a plot like this.

An example of corner detection with the Harris corner detector.

An overview of different approaches to corner detection, including improvements on the Harris detector and further developments, see e.g. http://en.wikipedia.org/wiki/Corner_detection.