Implementing a Harris corner detector

Question

I am implementing a Harris corner detector for educational purposes but I'm stuck at the harris response part. Basically, what I am doing, is:

1. Compute image intensity gradients in x- and y-direction
2. Blur output of (1)
3. Compute Harris response over output of (2)
4. Suppress non-maximas in output of (3) in a 3x3-neighborhood and threshold output

1 and 2 seem to work fine; however, I get very small values as the Harris response, and no point does reach the threshold. Input is a standard outdoor photography.

[...]
[Ix, Iy] = intensityGradients(img);
g = fspecial('gaussian');
Ix = imfilter(Ix, g);
Iy = imfilter(Iy, g);
H = harrisResponse(Ix, Iy);
[...]

function K = harrisResponse(Ix, Iy)
    max = 0;
    [sy, sx] = size(Ix);
    K = zeros(sy, sx);
    for i = 1:sx,
        for j = 1:sy,
            H = [Ix(j,i) * Ix(j,i), Ix(j,i) * Iy(j,i)
                Ix(j,i) * Iy(j,i), Iy(j,i) * Iy(j,i)];
            K(j,i) = det(H) / trace(H);
            if K(j,i) > max,
                max = K(j,i);
            end
        end
    end
    max
end

For the sample picture, max ends up being 6.4163e-018 which seems far too low.

Solution 1

A corner in Harris corner detection is defined as "the highest value pixel in a region" (usually 3X3 or 5x5) so your comment about no point reaching a "threshold" seems strange to me. Just collect all pixels that have a higher value than all other pixels in the 5x5 neighborhood around them.

Apart from that: I'm not 100% sure, but I think you should have:

K(j,i) = det(H) - lambda*(trace(H)^2) Where lambda is a positive constant that works in your case (and Harris suggested value is 0.04).

In general the only sensible moment to filter your input is before this point:

[Ix, Iy] = intensityGradients(img);

Filtering Ix2, Iy2 and Ixy doesn't make much sense to me.

Further, I think your sample code is wrong here (does function harrisResponse have two or three input variables?):

H = harrisResponse(Ix2, Ixy, Iy2);
[...]

function K = harrisResponse(Ix, Iy)

Solution 2

The solution that I implemented with python, it works for me I hope you find what you are looking for

import numpy as np
import matplotlib.pyplot as plt
from PIL.Image import *
from scipy import ndimage

def imap1(im):
    print('testing the picture . . .')
    a = Image.getpixel(im, (0, 0))
    if type(a) == int:
        return im
    else:
        c, l = im.size
        imarr = np.asarray(im)
        neim = np.zeros((l, c))
        for i in range(l):
            for j in range(c):
                t = imarr[i, j]
                ts = sum(t)/len(t)
                neim[i, j] = ts
        return neim

def Harris(im):
    neim = imap1(im)
    imarr = np.asarray(neim, dtype=np.float64)
    ix = ndimage.sobel(imarr, 0)
    iy = ndimage.sobel(imarr, 1)
    ix2 = ix * ix
    iy2 = iy * iy
    ixy = ix * iy
    ix2 = ndimage.gaussian_filter(ix2, sigma=2)
    iy2 = ndimage.gaussian_filter(iy2, sigma=2)
    ixy = ndimage.gaussian_filter(ixy, sigma=2)
    c, l = imarr.shape
    result = np.zeros((c, l))
    r = np.zeros((c, l))
    rmax = 0
    for i in range(c):
        print('loking for corner . . .')
        for j in range(l):
            print('test ',j)
            m = np.array([[ix2[i, j], ixy[i, j]], [ixy[i, j], iy2[i, j]]], dtype=np.float64)
            r[i, j] = np.linalg.det(m) - 0.04 * (np.power(np.trace(m), 2))
            if r[i, j] > rmax:
                rmax = r[i, j]
    for i in range(c - 1):
        print(". .")
        for j in range(l - 1):
            print('loking')
            if r[i, j] > 0.01 * rmax and r[i, j] > r[i-1, j-1] and r[i, j] > r[i-1, j+1]\
                                     and r[i, j] > r[i+1, j-1] and r[i, j] > r[i+1, j+1]:
                result[i, j] = 1

    pc, pr = np.where(result == 1)
    plt.plot(pr, pc, 'r+')
    plt.savefig('harris_test.png')
    plt.imshow(im, 'gray')
    plt.show()
    # plt.imsave('harris_test.png', im, 'gray')

im = open('chess.png')
Harris(im)

Solution 3

Basically, Harris corner detection will have 5 steps:

1. Gradient computation
2. Gaussian smoothing
3. Harris measure computation
4. Non-maximum suppression
5. Thresholding

If you are implementing in MATLAB, it will be easy to understand the algorithm and get the results.

The following code of MATLAB may help you to solve your doubts:

% Step 1: Compute derivatives of image
Ix = conv2(im, dx, 'same');
Iy = conv2(im, dy, 'same');

% Step 2: Smooth space image derivatives (gaussian filtering)
Ix2 = conv2(Ix .^ 2, g, 'same');
Iy2 = conv2(Iy .^ 2, g, 'same');
Ixy = conv2(Ix .* Iy, g, 'same');

% Step 3: Harris corner measure
harris = (Ix2 .* Iy2 - Ixy .^ 2) ./ (Ix2 + Iy2);

% Step 4: Find local maxima (non maximum suppression)
mx = ordfilt2(harris, size .^ 2, ones(size));

% Step 5: Thresholding
harris = (harris == mx) & (harris > threshold);

Solution 4

Proposed implementation is terribly inefficient. Lets' start after calculating gradients (which can be optimized too):

A = Ix.^2;
B = Iy.^2;
C = (Ix.*Iy).^4;
lambda = 0.04;

H = (A.*B - C) - lambda*(A+B).^2;

% if you really need max:
max(H(:))

No loops required, because Matlab hates loops.

View more solutions

Comments

[Etan]

I have reverted to not filter Ix2 etc anymore, therefore there was some bug left in the copy on stackoverflow.

[Etan]

The problem was that I didn't sum up all the pixels in the 3x3 square to find out the Ix2 etc; instead, I have just used the corresponding pixel. After changing H in a way that it sums up all Ix2, Ixy and Iy2 for all 9 pixels, it looks very nice.

[Etan]

det(H)/trace(H) is a ofen-used approximation in the case where you won't have a lambda.

[jilles de wit]

I didn't know about that last trick. Nice. It seems you solved the problem yourself, nice! (and the age old trick still works: just explaining the problem to someone helps you get to the solution) is this the working code ?

[Hilder Vitor Lima Pereira]

But why calculate C = (Ix.*Iy).^4 instead of simple C = (Ix.*Iy) ?

[acoustic python]

what is pc, and pr what it gives?

[Walid Bousseta]

@acousticpython pc and pr are index where result == 1, means that result[pc][pr] == 1, items in pc and pr at the same position are where result equal to one in the 2d array