Experimenting with the Harris Corner Detector Algorithm in Matlab

tl;dr

Example code to start experimenting with the Harris Corner Detection Algorithm on your own, in Matlab. And some of the results I obtained in my own testing.

Corner Cases

Among the classic algorithms in Computer Vision is Harris Corner Detection. (For the original paper, from 1988, see here.) The problem it solves: Given any image file, can you find which areas correspond to a corner with a high degree of certainty?

The answer is yes. And the algorithm is elegant. Essentially the steps are:

Convert your image to a Matrix, storing the pixel values by intensity (converting to B/W makes this easier)
Find the horizontal and vertical gradients of the image (essentially the derivative of pixel intensities). This will let us know where there are horizontal and vertical edges.
Scan the image getting groups of pixels for however wide an area you are concerned with, and perform a calculation to determine where the horizontal and vertical edges intersect. (This requires some relatively complex linear algebra, but I’m not focusing on the theory here, just the implementation.)

Recently, I was trying to implement my own version of the Harris Detector in Matlab, and ended up banging my head against the wall for a few hours while figuring out some of the subtler details the hard way.

Here’s the code I came up with, and some examples of the outputs. Grab the code from both gists below, and you can start experimenting on your own.

Harris Algorithm

	function [ x, y, scores, Ix, Iy ] = harris_corners( image )
	%HARRIS_CORNERS Extracts points with a high degree of 'cornerness' from
	%RGB image matrix of type uint8
	% Input – image = NxMx3 RGB image matrix
	% Output – x = nx1 vector denoting the x location of each of n
	% detected keypoints
	% y = nx1 vector denoting the y location of each of n
	% detected keypoints
	% scores = an nx1 vector that contains the value (R) to which a
	% a threshold was applied, for each keypoint
	% Ix = A matrix with the same number of rows and columns as the
	% input image, storing the gradients in the x-direction at each
	% pixel
	% Iy = A matrix with the same nuimber of rwos and columns as the
	% input image, storing the gradients in the y-direction at each
	% pixel

	% compute the gradients, re-use code from HW2P, use window size of 5px
	% convert image to grayscale first
	G = rgb2gray(image);

	% convert to double
	G2 = im2double(G);

	% create X and Y Sobel filters
	horizontal_filter = [1 0 –1; 2 0 –2; 1 0 –1];
	vertical_filter = [1 2 1; 0 0 0 ; –1 –2 –1];

	% using imfilter to get our gradient in each direction
	filtered_x = imfilter(G2, horizontal_filter);
	filtered_y = imfilter(G2, vertical_filter);

	% store the values in our output variables, for clarity
	Ix = filtered_x;
	Iy = filtered_y;

	% Compute the values we need for the matrix…
	% Using a gaussian blur, because I get more positive values after applying
	% it, my values all skew negative for some reason…
	f = fspecial('gaussian');
	Ix2 = imfilter(Ix.^2, f);
	Iy2 = imfilter(Iy.^2, f);
	Ixy = imfilter(Ix.*Iy, f);

	% set empirical constant between 0.04-0.06
	k = 0.04;

	num_rows = size(image,1);
	num_cols = size(image,2);

	% create a matrix to hold the Harris values
	H = zeros(num_rows, num_cols);

	% % get our matrix M for each pixel
	for y = 6:size(image,1)-6 % avoid edges
	for x = 6:size(image,2)-6 % avoid edges
	% calculate means (because mean is sum/num pixels)
	% generally, this algorithm calls for just finding a sum,
	% but using the mean makes visualization easier, in my code,
	% and it doesn't change which points are computed to be corners.
	% Ix2 mean
	Ix2_matrix = Ix2(y–2:y+2,x–2:x+2);
	Ix2_mean = sum(Ix2_matrix(:));

	% Iy2 mean
	Iy2_matrix = Iy2(y–2:y+2,x–2:x+2);
	Iy2_mean = sum(Iy2_matrix(:));

	% Ixy mean
	Ixy_matrix = Ixy(y–2:y+2,x–2:x+2);
	Ixy_mean = sum(Ixy_matrix(:));

	% compute R, using te matrix we just created
	Matrix = [Ix2_mean, Ixy_mean;
	Ixy_mean, Iy2_mean];
	R1 = det(Matrix) – (k * trace(Matrix)^2);

	% store the R values in our Harris Matrix
	H(y,x) = R1;

	end
	end

	% set threshold of 'cornerness' to 5 times average R score
	avg_r = mean(mean(H));
	threshold = abs(5 * avg_r);

	[row, col] = find(H > threshold);

	scores = [];
	%get all the values
	for index = 1:size(row,1)
	%see what the values are
	r = row(index);
	c = col(index);
	scores = cat(2, scores,H(r,c));
	end

	y = row;
	x = col;

	end

view raw

harris_corners.m

hosted with ❤ by GitHub

Displaying Your Results

	image = imread('your_image.png');
	image = imresize(image, 0.75); # or however much you want to resize if your image is large
	[ x, y, scores, Ix, Iy ] = harris_corners( image );
	figure; imshow(image)
	hold on
	for i = 1:size(scores,2)
	plot(x(i), y(i), 'ro', 'MarkerSize', scores(i) * 2); # you may need to play with this multiplier or divisor based on your image
	# I've used –> (/1000) to (* 10)

	end
	saveas(gcf,'your_image_with_corners.png');
	hold off

view raw

display_corners.m

hosted with ❤ by GitHub

My Outputs

Here’s what you can expect to see. I’m coloring the corners RED, using circles that grow larger the more likely the area is to hold a corner. So, larger circles –> more obvious corners. (Also note that I’m downsizing my output images slightly, just to make computation faster. )

Closing

That’s it. With this code, you can start experimenting with Harris Detectors in about 30 seconds. And hopefully the way I’ve written it is transparent enough to make it obvious what’s going on behind the scenes. Most of the implementations out there are somewhat opaque and hard to understand at a conceptual level. But I’m hoping that my contribution will make this all a little more straightforward to those just gettings started. Enjoy. -t

tl;dr

Corner Cases

Harris Algorithm

Displaying Your Results

My Outputs

Closing

Share this:

Like this:

Related

Published by adpoe