classdef FeatureDetector
% General wrapper for different feature detectors. The options of the
% The parameters can be changes by accessing the Opts structure.
% The following methods (detectors) are available:
% - Harris
% - Hessian
% - Hessian Squared - (Hessian2)
% - SURF
% - BRISK
% - SIFT (requires VL_Fleat SIFT installation)
%
% =========================================================================
%
% Copyright (c) 2016 Anders Hast, Damian Matuszewski
% Centre for Image Analysis
% Uppsala University
% 
% 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:
%
%  1) Any scientific publications coming from using or modifying this code
%     should cite the original paper:
%       Matuszewski, D.J., Hast, A., Wählby, C., Sintorn, I.-M. (2016)
%       A short feature vector for image matching: the Log Polar Magnitude 
%       feature descriptor. Computer Vision and Image Understanding [submitted]       
%
%  2) The above copyright notice and this permission notice shall be 
%     included in all copies or substantial portions of the Software.
%
%  3) The Software is provided "as is", without warranty of any kind.
%
% =========================================================================
%
  properties (SetAccess=public, GetAccess=public)
    % The default values of the detectors options
    Opts = struct(...
            'Method', 'SURF', ... % Name of the detector     
            'VL_Feat_Path', '', ... % Path to the .../vlfeat-X.X.XX/toolbox directory in the VL_Feat SIFT installation 
            'PointsNo', 1000, ... % Number of feature points (only for Hessian and Harris)
            'MinDistBetweenPoints', 1, ... % (only for Hessian and Harris)
            'FilterKernel', 'Cubic', ... % (only for Hessian and Harris)
            'PreSmoothGausSigma', 2.0 ... % (only for Hessian and Harris)
      );
  end
  
  methods
    function obj = FeatureDetector()
      % Constructor: add source code to the MATLAB path.

      % Setup the paths
      addpath(fullfile('Detectors_code'));   
    end

    function [frames] = getFrames(obj, imagePath)
      rads = [];
        
      % Get the input image
      I = imread(imagePath);
      info = imfinfo(imagePath);
      imgSize = size(I);
      imgChannels = 1;
      if length(imgSize) > 2
        if imgSize(3) > 1
            I = rgb2gray(I);
            imgChannels = imgSize(3);
        end
      end
      im = double(I)/(2^(info.BitDepth/imgChannels) - 1);
        
      % presmooth the image
      hx = fspecial('gaussian',[7,1],1);
      hy = fspecial('gaussian',[1,7],1);
      im = conv2(hx,hy,im,'same');  % blur image with sigma=1.0

      % initialize the filter krernels for Hessian and Harris
      [k, d, d2]=kernel(obj.Opts.FilterKernel);    
      g2 = fspecial('gaussian', obj.Opts.PreSmoothGausSigma*6 + 1, obj.Opts.PreSmoothGausSigma); % Gaussian with sigma  
      
      % get the frames with chosen detector
      if strcmp(obj.Opts.Method, 'Harris') % Harris (corners)
        Iy  =-conv2(conv(d,k),conv(k,k),im,'same');
        Ix  =-conv2(conv(k,k),conv(d,k),im,'same');
        Ix2  = conv2(Ix.*Ix,g2,'same');
        Iy2  = conv2(Iy.*Iy,g2,'same');
        IxIy = conv2(Ix.*Iy,g2,'same');

        cim = (Ix2.*Iy2 - IxIy.^2)./(Ix2 + Iy2 + eps); % My preferred  measure.
        [rows, cols] = nonmaxsupp(cim, obj.Opts.PointsNo, obj.Opts.MinDistBetweenPoints, 0);
      elseif strcmp(obj.Opts.Method, 'Hessian') % Hessian (blobs)
        Ixx = conv2(conv(k,k),conv(d2,k),im,'same');
        Iyy = conv2(conv(d2,k),conv(k,k),im,'same'); %
        Ixy = conv2(conv(d,k),conv(d,k),im,'same');

        Ixx = conv2(Ixx,g2,'same');
        Iyy = conv2(Iyy,g2,'same');
        Ixy = conv2(Ixy,g2,'same');

        cim = (Ixx.*Iyy - Ixy.^2);
        [rows, cols] = nonmaxsupp(cim, obj.Opts.PointsNo, obj.Opts.MinDistBetweenPoints, 0);            
      elseif strcmp(obj.Opts.Method, 'Hessian2') % Hessian squared (blobs)
        Ixx = conv2(conv(k,k),conv(d2,k),im,'same');
        Iyy = conv2(conv(d2,k),conv(k,k),im,'same'); %
        Ixy = conv2(conv(d,k),conv(d,k),im,'same');

        Ixx = conv2(Ixx,g2,'same');
        Iyy = conv2(Iyy,g2,'same');
        Ixy = conv2(Ixy,g2,'same');

        cim = (Ixx.*Iyy - Ixy.^2).^2;
        [rows, cols] = nonmaxsupp(cim, obj.Opts.PointsNo, obj.Opts.MinDistBetweenPoints, 0); 
      elseif strcmp(obj.Opts.Method, 'SURF') % SURF features (blobs)
        points = detectSURFFeatures(im);
        rows = round(points.Location(:,2));
        cols = round(points.Location(:,1));
        rads = points.Scale;
      elseif strcmp(obj.Opts.Method, 'BRISK') % BRISK features (corners)
        points = detectBRISKFeatures(im);
        rows = round(points.Location(:,2));
        cols = round(points.Location(:,1));
        rads = points.Scale;
      elseif strcmp(obj.Opts.Method, 'SIFT') % SIFT features (corners) requires VL_feat installation
        run([obj.Opts.VL_Feat_Path, '/vlfeat-0.9.20/toolbox/vl_setup']) % initialize VL_Feat SIFT
        I = single(I);
        p = vl_sift(I); %,'PeakThresh', 4);        
        rows = p(2,:)';
        cols = p(1,:)';
        rads = p(3,:)';
      else
        error('\nWrong feature detector name! \nPlease, use one of the following: Harris, Hessian, Hessian2, SURF, BRISK, or SIFT (requires VL_Feat SIFT installed).')
      end
      fprintf('\n%d features in %s detected using %s.', length(cols), imagePath, obj.Opts.Method) 
      
      if isempty(rads)
        frames = [cols, rows]';
      else
        frames = [cols, rows, rads]';
      end
      
    end
    
  end  
end