% 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 nonmaxsupptsfunction [cim, r, c, rsubp, csubp] = harris(im, sigma, thresh, radius, disp)        error(nargchk(2,5,nargin));    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