Corner Detection using Harris

Harris and Stephens improved upon Moravec's corner detector by considering the differential of the corner score with respect to direction directly, instead of using shifted patches. (This corner score is often referred to as autocorrelation, since the term is used in the paper in which this detector is described. However, the mathematics in the paper clearly indicate that the sum of squared differences is used.)

C. Harris and M. Stephens (1988). "A combined corner and edge detector" (PDF). Proceedings of the 4th Alvey Vision Conference. pp. 147--151. http://www.csse.uwa.edu.au/~pk/research/matlabfns/Spatial/Docs/Harris/A_Combined_Corner_and_Edge_Detector.pdf.

How does it work?
The first step in Harris is to compute a corner response function. Harris uses a series of filters that are considered the most important steps to find a corner. These kernels are:

• A presmoothing filter pfilt = {0.223755f,0.552490f,0.223755f};
• a gradient filter gfilt = {0.453014f,0.0f,-0.453014f};
• a Blur Filter bfilt = {0.01563f,0.09375f,0.234375f,0.3125f, 4 0.234375f,0.09375f,0.01563f};

The combination of these filters and few parameters are what make Harris successful in finding corners. The parameters are as followed::

• Steering parameter 0.04 - 0.06 (0.05 default)
• Response threshold 10K - 1M (default 25,000)
• neighborhood radius 10 pixels

clear all
close all
clc;

disp('Harris Corner detection Purposes Only....');

fin = 'smallVersion2.avi'; %input File name to be changed
fout = 'movie2.avi'; %output File name to be changed

%get the File info
fileinfo = aviinfo(fin);
%get the number of Frames
nframes = fileinfo.NumFrames;

aviobj = avifile(fout, 'compression', 'none', 'fps',fileinfo.FramesPerSecond);

for i = 1:nframes %You may need to limit the nFrames, may result in a large file size
%Read frames from input video
mov_in = aviread(fin,i);
im_in = frame2im(mov_in);

%Do processing on each frame of the video

im = rgb2gray(im_in);
% Find Harris corners in image1 and image2
if ~isa(im,'double')
im = double(im);
end

subpixel = nargout == 5;

dx = [-1 0 1; -1 0 1; -1 0 1]; % Derivative masks
dy = dx';

Ix = conv2(im, dx, 'same'); % Image derivatives
Iy = conv2(im, dy, 'same');

% Generate Gaussian filter of size 6*sigma (+/- 3sigma) and of
% minimum size 1x1.
sigma =1;
g = fspecial('gaussian',max(1,fix(6*sigma)), sigma);

Ix2 = conv2(Ix.^2, g, 'same'); % Smoothed squared image derivatives
Iy2 = conv2(Iy.^2, g, 'same');
Ixy = conv2(Ix.*Iy, g, 'same');
k = 0.04;
cim = (Ix2.*Iy2 - Ixy.^2) - k*(Ix2 + Iy2).^2; % Original Harris measure.


if nargin > 2 % We should perform nonmaximal suppression and threshold

if subpixel
[r,c,rsubp,csubp] = nonmaxsuppts(cim, 1, 10000, im);
else
[r,c] = nonmaxsuppts(cim, 1, 10000, im);
end
end


show(im,1), hold on, plot(c,r,'r+');
drawnow


%Write frames to output video
%frm = im2frame(im_out);
%aviobj = addframe(aviobj,frm);

%i %Just to display frame number you can ommit
end;

%Don't forget to close output file
aviobj = close(aviobj);
return;




function [r,c, rsubp, csubp] = nonmaxsuppts(cim, radius, thresh, im)

subPixel = nargout == 4; % We want sub-pixel locations
[rows,cols] = size(cim);

% Extract local maxima by performing a grey scale morphological
% dilation and then finding points in the corner strength image that
% match the dilated image and are also greater than the threshold.

sze = 2*radius+1; % Size of dilation mask.
mx = ordfilt2(cim,sze^2,ones(sze)); % Grey-scale dilate.

% Make mask to exclude points within radius of the image boundary.
bordermask = zeros(size(cim));
bordermask(radius+1:end-radius, radius+1:end-radius) = 1;

% Find maxima, threshold, and apply bordermask
cimmx = (cim==mx) & (cim>thresh) & bordermask;

[r,c] = find(cimmx); % Find row,col coords.


if subPixel % Compute local maxima to sub pixel accuracy
if ~isempty(r) % ...if we have some ponts to work with

ind = sub2ind(size(cim),r,c); % 1D indices of feature points
w = 1; % Width that we look out on each side of the feature
% point to fit a local parabola

% Indices of points above, below, left and right of feature point
indrminus1 = max(ind-w,1);
indrplus1 = min(ind+w,rows*cols);
indcminus1 = max(ind-w*rows,1);
indcplus1 = min(ind+w*rows,rows*cols);

% Solve for quadratic down rows
cy = cim(ind);
ay = (cim(indrminus1) + cim(indrplus1))/2 - cy;
by = ay + cy - cim(indrminus1);
rowshift = -w*by./(2*ay); % Maxima of quadradic

% Solve for quadratic across columns
cx = cim(ind);
ax = (cim(indcminus1) + cim(indcplus1))/2 - cx;
bx = ax + cx - cim(indcminus1);
colshift = -w*bx./(2*ax); % Maxima of quadradic

rsubp = r+rowshift; % Add subpixel corrections to original row
csubp = c+colshift; % and column coords.
else
rsubp = []; csubp = [];
end
end

if nargin==4 & ~isempty(r) % Overlay corners on supplied image.
figure(1), imshow(im,[]), hold on
if subPixel
plot(csubp,rsubp,'r+'), title('corners detected');
else
plot(c,r,'r+'), title('corners detected');
end
end