harris.m

用matlab实现Harris Corner Detector的算法

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Final year Project Simulation Code                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Navigation Parameter Estimation of Vision Guided Helicopter           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Members       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Kashif Shahzad    01-ET-31
% Imran Raees       01-ET-76
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Harris interset point detector %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------------------------------------------------------------------------%
% This Program calculates the Feature points in an image/frame based
% on the research by HARRIS
% Syntax [cornerness, r, c] = harris(im, sigma, thresh, radius, disp)
%-------------------------------------------------------------------------%

function [cornerness, r, c] = harris(im, sigma, thresh, radius, disp)
% cornerness = Corner values
% r = rows (x coordinate of feature point)
% c = columns (y coordinate of feature points)

% number of arguments check and error msg
error(nargchk(2,5,nargin)); 
    
% Derivative masks
dy = [-1 0 1; -1 0 1; -1 0 1]; 
dx = dy'; %dx is the transpose matrix of dy

% Ix and Iy are the horizontal and vertical edges of image
Ix = conv2(im, dx, 'same');    
Iy = conv2(im, dy, 'same');    

% Calculating the gradient of the image Ix and Iy
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');

% My preferred  measure according to research paper
cornerness = (Ix2.*Iy2 - Ixy.^2)./(Ix2 + Iy2 + eps); 

% We should perform nonmaximal suppression and threshold
if nargin > 2 
    sze = 2*radius+1;                                  % Size of mask
	mx = ordfilt2(cornerness,sze^2,ones(sze));         % Grey-scale dilate
	cornerness = (cornerness==mx)&(cornerness>thresh); % Find maxima
	[r,c] = find(cornerness);                          % Find row,col coords.
else
    r = []; c = [];
end