%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