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.

Answer 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)

Answer 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)