function points = kp_harris(im)
% Extract keypoints using Harris algorithm (with an improvement
% version)
% INPUT
% =====
% im : the graylevel image
%
% OUTPUT
% ======
% points : the interest points extracted
%
% REFERENCES
% ==========
% C.G. Harris and M.J. Stephens. "A combined corner and edge detector",
% Proceedings Fourth Alvey Vision Conference, Manchester.
% pp 147-151, 1988.
%
% Alison Noble, "Descriptions of Image Surfaces", PhD thesis, Department
% of Engineering Science, Oxford University 1989, p45.
%
% C. Schmid, R. Mohrand and C. Bauckhage, "d",
% Int. Journal of Computer Vision, 37(2), 151-172, 2000.
%
% EXAMPLE
% =======
% points = kp_harris(im)
% only luminance value
%size(im)
im = double(im(:,:,1));
sigma = 1.5;
% derivative masks
s_D = 0.7*sigma;
x = -round(3*s_D):round(3*s_D);
dx = x .* exp(-x.*x/(2*s_D*s_D)) ./ (s_D*s_D*s_D*sqrt(2*pi));
dy = dx';
% image derivatives
Ix = conv2(im, dx, 'same');
Iy = conv2(im, dy, 'same');
% sum of the Auto-correlation matrix
s_I = sigma;
g = fspecial('gaussian',max(1,fix(6*s_I+1)), s_I);
Ix2 = conv2(Ix.^2, g, 'same'); % Smoothed squared image derivatives
Iy2 = conv2(Iy.^2, g, 'same');
Ixy = conv2(Ix.*Iy, g, 'same');
% interest point response
cim = (Ix2.*Iy2 - Ixy.^2)./(Ix2 + Iy2 + eps);
% find local maxima on 3x3 neighborgood
[r,c,max_local] = findLocalMaximum(cim,3*s_I);
% set threshold 1% of the maximum value
%t = 0.01*max(max_local(:));
t = 0.6*max(max_local(:)); %door.jpg
%t = 0.48*max(max_local(:)); %sunflower.jpg
% find local maxima greater than threshold
[r,c] = find(max_local>=t);
% build interest points
points = [r,c];
end
function [row,col,max_local] = findLocalMaximum(val,radius)
% Determine the local maximum of a given value
%
%
% INPUT
% =====
% val : the NxM matrix containing values
% radius : the radius of the neighborhood
%
% OUTPUT
% ======
% row : the row position of the local maxima
% col : the column position of the local maxima
% max_local : the NxM matrix containing values of val on unique local maximum
%
% EXAMPLE
% =======
% [l,c,m] = findLocalMaximum(img,radius);
% FIND UNIQUE LOCAL MAXIMA USING FILTERING (FAST)
mask = fspecial('disk',radius)>0;
nb = sum(mask(:));
highest = ordfilt2(val, nb, mask);
second_highest = ordfilt2(val, nb-1, mask);
index = highest==val & highest~=second_highest;
max_local = zeros(size(val));
max_local(index) = val(index);
[row,col] = find(index==1);
end
结果: