%IM_HARRIS Fixed mapping executing the Harris corner detector
%
% X = IM_HARRIS(A,N,SIGMA)
% X = A*IM_HARRIS([],N,SIGMA)
% X = A*IM_HARRIS(N,SIGMA)
%
% INPUT
% A Datafile or dataset with images
% N Number of desired Harris points per image (default 100)
% SIGMA Smoothing size (default 3)
%
% OUTPUT
% X Dataset with a [N,3] array with for every image
% x, y and strength per Harris point.
%
% DESCRIPTION
% We use Kosevi's [1] software to find the corner points according to
% Harris [2]. On top of the Kosevi Harris point detector we run
% - multi-feature images (e.g. color images) are averaged
% - only points that are maximum in a K x K window are selected.
% If less points are found, K is iteratively reduced.
% The initial value of K is about 4*SIGMA.
% Although SIGMA can be interpreted as scaling parameter, it might be
% better to appropriately subsample images instead of using a large SIGMA.
%
% If you use this software for publications, please refer to [1] and [2].
%
% REFERENCES
% [1] P. D. Kovesi, MATLAB and Octave Functions for Computer Vision and
% Image Processing, School of Computer Science & Software Engineering,
% The University of Western Australia. Available from:
% <http://www.csse.uwa.edu.au/~pk/research/matlabfns/>.
% [2] C. Harris and M. Stephens, A combined corner and edge detector,
% Proc. 4th Alvey Vision Conf., 1988, pp. 147-151.
%
% EXAMPLE
% delfigs
% a = kimia; % take simple shapesas example
% b = gendat(a,25)*im_gray; % just 25 images at random
% c = data2im(b); % convert dataset to images for display
% x = im_harris(b,15,1); % compute maximum 15 Harris points at scale 1
% y = data2im(x); % unpack dataset with results
% for j=1:25 % show results one by one
% figure(j); imagesc(c(:,:,1,j)); colormap gray; hold on
% scatter(y(:,1,1,j),y(:,2,1,j),'r*');
% end
% showfigs
function x = im_harris(varargin)
argin = shiftargin(varargin,'scalar');
argin = setdefaults(argin,[],100,3,[]);
if mapping_task(argin,'definition')
x = define_mapping(argin,'fixed');
x = setname(x,'Harris points');
else
[a,n,s,k] = deal(argin{:});
if isa(a,'prdataset') % allows datafiles too
isobjim(a);
a = im_gray(a);
x = filtim(a,mfilename,{n,s});
else
a = double(a);
if size(a,3) ~= 1
error('2D gray value image expected');
end
[mr,mc] = size(a); % here we put a border of 15 pixels around the image
aa = bord(a,NaN,15); % mirroring the border to avoid artificial corners
cc = harris(aa,s); % at the image border, finally we take the
c = resize(cc,15,mr,mc); % original image size (bord and resize are in ./private)
fsizes = [25,21,17,13,11,9,7,5,3];
F = find(fsizes < 4*s);
fsizes = fsizes(F);
for fsize = fsizes
shape = ones(fsize,fsize);
cmax = ordfilt2(c,fsize*fsize,shape);
J = find(cmax == c); % are a local maximum
[Iy,Ix] = ind2sub(size(c),J); % second attempt to get rid of image border points
K = find(Iy > 3 & Ix >3 & Iy < (mr-2) & Ix < (mc-2));
J = J(K);
Z = find(cmax(J) > 0); % find non-zeros only
J = J(Z);
if length(J) > n % if we have sufficient points
break; % stop
end % otherwise, get more local maxima
end
cs = cmax(J);
[ss,R] = sort(-cs); % rank Harris points
if length(J) >= n % find x,y coordinates of first n Harris points
[Iy,Ix] = ind2sub(size(c),J(R(1:n)));
x = [Ix, Iy, -ss(1:n)];
else
x = zeros(n,3);
[Iy,Ix] = ind2sub(size(c),J(R));
x(1:length(R),:) = [Ix, Iy, -ss(1:length(R))];
end
end
end
return
% HARRIS - Harris corner detector
%
% Usage: cim = harris(im, sigma)
% [cim, r, c] = harris(im, sigma, thresh, radius, disp)
% [cim, r, c, rsubp, csubp] = harris(im, sigma, thresh, radius, disp)
%
% Arguments:
% im - image to be processed.
% sigma - standard deviation of smoothing Gaussian. Typical
% values to use might be 1-3.
% thresh - threshold (optional). Try a value ~1000.
% radius - radius of region considered in non-maximal
% suppression (optional). Typical values to use might
% be 1-3.
% disp - optional flag (0 or 1) indicating whether you want
% to display corners overlayed on the original
% image. This can be useful for parameter tuning. This
% defaults to 0
%
% Returns:
% cim - binary image marking corners.
% r - row coordinates of corner points.
% c - column coordinates of corner points.
% rsubp - If five return values are requested sub-pixel
% csubp - localization of feature points is attempted and
% returned as an additional set of floating point
% coords. Note that you may still want to use the integer
% valued coords to specify centres of correlation windows
% for feature matching.
%
% If thresh and radius are omitted from the argument list only 'cim' is returned
% as a raw corner strength image. You may then want to look at the values
% within 'cim' to determine the appropriate threshold value to use. Note that
% the Harris corner strength varies with the intensity gradient raised to the
% 4th power. Small changes in input image contrast result in huge changes in
% the appropriate threshold.
% 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.
% Copyright (c) 2002-2005 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.
% March 2002 - original version
% December 2002 - updated comments
% August 2005 - changed so that code calls nonmaxsuppts
function [cim, r, c, rsubp, csubp] = harris(im, sigma, thresh, radius, disp)
narginchk(2,5);
if nargin == 4
disp = 0;
end
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.
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');
% Compute the Harris corner measure. Note that there are two measures
% that can be calculated. I prefer the first one below as given by
% Nobel in her thesis (reference above). The second one (commented out)
% requires setting a parameter, it is commonly suggested that k=0.04 - I
% find this a bit arbitrary and unsatisfactory.
cim = (Ix2.*Iy2 - Ixy.^2)./(Ix2 + Iy2 + eps); % My preferred measure.
% 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 disp % Call nonmaxsuppts to so that image is displayed
if subpixel
[r,c,rsubp,csubp] = nonmaxsuppts(cim, radius, thresh, im);
else
[r,c] = nonmaxsuppts(cim, radius, thresh, im);
end
else % Just do the nonmaximal suppression
if subpixel
[r,c,rsubp,csubp] = nonmaxsuppts(cim, radius, thresh);
else
[r,c] = nonmaxsuppts(cim, radius, thresh);
end
end
end
return
% NONMAXSUPPTS - Non-maximal suppression for features/corners
%
% Non maxima suppression and thresholding for points generated by a feature
% or corner detector.
%
% Usage: [r,c] = nonmaxsuppts(cim, radius, thresh, im)
% /
% optional
%
% [r,c, rsubp, csubp] = nonmaxsuppts(cim, radius, thresh, im)
%
% Arguments:
% cim - corner strength image.
% radius - radius of region considered in non-maximal
% suppression. Typical values to use might
% be 1-3 pixels.
% thresh - threshold.
% im - optional image data. If this is supplied the
% thresholded corners are overlayed on this
% image. This can be useful for parameter tuning.
% Returns:
% r - row coordinates of corner points (integer valued).
% c - column coordinates of corner points.
% rsubp - If four return values are requested sub-pixel
% csubp - localization of feature points is attempted and
% returned as an additional set of floating point
% coords. Note that you may still want to use the integer
% valued coords to specify centres of correlation windows
% for feature matching.
%
% Note: An issue with integer valued images is that if there are multiple pixels
% all with the same value within distance 2*radius of each other then they will
% all be marked as local maxima.
% Copyright (c) 2003-2013 Peter Kovesi
% Centre for Exploration Targeting
% The University of Western Australia
% peter.kovesi at uwa edu au
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.
% September 2003 Original version
% August 2005 Subpixel localization and Octave compatibility
% January 2010 Fix for completely horizontal and vertical lines (by Thomas Stehle,
% RWTH Aachen University)
% January 2011 Warning given if no maxima found
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
rowshift = zeros(size(ind));
cy = cim(ind);
ay = (cim(indrminus1) + cim(indrplus1))/2 - cy;
by = ay + cy - cim(indrminus1);
rowshift(ay ~= 0) = -w*by(ay ~= 0)./(2*ay(ay ~= 0)); % Maxima of quadradic
rowshift(ay == 0) = 0;
% Solve for quadratic across columns
colshift = zeros(size(ind));
cx = cim(ind);
ax = (cim(indcminus1) + cim(indcplus1))/2 - cx;
bx = ax + cx - cim(indcminus1);
colshift(ax ~= 0) = -w*bx(ax ~= 0)./(2*ax(ax ~= 0)); % Maxima of quadradic
colshift(ax == 0) = 0;
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
hold off
end
if isempty(r)
warning('No maxima above threshold found\n');
end
return