classdef LPMDescriptor
% LPMDescriptor - Log Polar Magnitude frame descriptor.

%   LPMDescriptor Constructs new Log Polar Magnitude feature descriptor. 
%   This algorithm requires precomputed frames (regions) in the form of discs:
%       frames(1,:) - x coordinate of the region center
%       frames(2,:) - y coordinate of the region center
%       frames(3,:) - (optional) scale or radius of the region 
%   
%   The algorithm first computes the Log Polar Transform in the region with
%   sampling resolution (number of rings and points per ring) defined in 
%   'Sampling' option. Then the log polar results are transformed to 
%   frequency domain with Fast Fourier Transform. Next, magnitudes of 
%   selected frequencies (defined as list of indexes in 'FreqSamplMaskIdx')
%   are used to build the descriptor. 
%
% Options:
%   FixedRadius:: false
%       Use fixed radius for all frames
%
%   Radius:: 32
%       Specifies the radius to be used in all frames
%
%   Sampling:: 32
%       Sampling resolution = number of rings and points per ring. Currently
%       only two values are supported: 16 | 32.
%
%   FrameScalingFactor:: 14
%       Scale convertion rate for handling different feature detectors. Use 
%       1 if the frames contain radius information in pixels. The scale /  
%       radius is multiplied by this scalar scaling factor before descriptor  
%       is calculated. This parameter is dataset specific! Modify it to get
%       better results for your data.
%       Some experimental values: SURF -> 9, SIFT -> 14, BRISK -> 1
%
%   FreqSamplMaskIdx:: []
%       Custom frequency mask to be used by descriptor. If not specified
%       the predefined masks for 16x16 or 32x32 sampling will be used
%       instead. This parameter is dataset specific! Modify it to get
%       better results for your data.
%
%
% =========================================================================
%
% 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 descriptor options
    Opts = struct(...
      'FixedRadius', false,... % use fixed radius for all frames
      'Radius', 32,... % the radius to be used in all frames
      'Sampling', 32,... % sampling resolution = number of rings and points per ring (16 | 32)
      'FrameScalingFactor', 14,... % scale convertion rate for different feature detectors (SURF -> 9, SIFT -> 14, BRISK -> 1)
      'FreqSamplMaskIdx', []... % custom frequency mask to be used by descriptor
      );
  end

  properties (Constant)
    % predefined frequency masks for building the descriptor with 32 x 32
    % sampling
    PredefFreqMask_16 = [17, 33, 49, 65, 81, 97, 18, 34, 50, 66, 82, 98,...
                       19, 35, 51, 67, 83, 99, 20, 36, 52, 68, 84, 100,...
                       145, 161, 177, 193, 209, 225, 146, 162, 178, 194,...
                       210, 226, 147, 163, 179, 195, 211, 227, 148, 164,...
                       180, 196, 212, 228];
    % predefined frequency masks for building the descriptor with 16 x 16
    % sampling
    PredefFreqMask_32 = [33, 65, 97, 129, 161, 193, 225, 34, 66, 98, 130,...
                        162, 194, 226, 35, 67, 99, 131, 163, 195, 227,...
                        36, 68, 100, 132, 164, 196, 228, 769, 801, 833,...
                        865, 897, 929, 961, 770, 802, 834, 866, 898, 930,...
                        962, 771, 803, 835, 867, 899, 931, 963, 772, 804,...
                        836, 868, 900, 932, 964];
    % image interpolation mask (Gaussian) for sampling
    SamplMaskSize = 5;
    % image interpolation sigma (Gaussian) for sampling
    SamplSigma = 1.0;                            
  end
  
  methods
    function obj = LPMDescriptor(varargin)
      % Constructor: covert .cpp source code to mex if necessary and add
      % source code to the MATLAB path.

      % Setup the paths
      addpath(fullfile('LPM_code'));   
      
      % Install mex files if necessary
      if exist('extractLPMFeatures', 'file') ~= 3
        mex LPM_code/extractLPMFeatures.cpp
      end
    end

    function [frames, descriptors] = extractDescriptors(obj, imagePath, frames)
      % EXTRACTDESCRIPTORS Compute Log Polar Magnitute feature vectors for disk frames.

      % set the frequency mask for building the descriptor
      if obj.Opts.Sampling == 16
        if isempty(obj.Opts.FreqSamplMaskIdx) % use predefined mask
          obj.Opts.FreqSamplMaskIdx = obj.PredefFreqMask_16;
        end
      elseif obj.Opts.Sampling == 32
        if isempty(obj.Opts.FreqSamplMaskIdx) % use predefined mask
          obj.Opts.FreqSamplMaskIdx = obj.PredefFreqMask_32;
        end 
      else
        error('Unsuported sampling resolution! Use ''Sampling'' = 16 or 32.');
      end
      
      % 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);
      
      % check if all frames are inside the allowed image area and fix their
      % scales
      [y, x, r] = obj.checkFeaturePoints(frames, imgSize); 
      frames = [x; y; r];

      % blur image with sigma=1.0
      hx = fspecial('gaussian',[7,1],1);
      hy = fspecial('gaussian',[1,7],1);
      im = conv2(hx,hy,im,'same'); 

      % calculate the feature vectors
      descriptors = extractLPMFeatures(im, y, x, r, obj.Opts.Sampling,...
                                      obj.SamplMaskSize, obj.SamplSigma,...
                                      double(obj.Opts.FreqSamplMaskIdx));          

      % normalize feature vectors to unit vectors
      for i = 1:length(descriptors)
          descriptors(:,i) = descriptors(:,i)./norm(descriptors(:,i),2);
      end      
      
      fprintf('\n%d feature vectors in %s described using LPM.', length(descriptors), imagePath)
    end   
  end

  methods(Access=protected)
    function [y, x, r] = checkFeaturePoints(obj, frames, imSize)
        % fixes the scaling and keeps only the frames that are entirely 
        % within the image      
        y = double(frames(2,:));
        x = double(frames(1,:));
        rows = imSize(1);
        cols = imSize(2);
        b = obj.SamplMaskSize;
        if obj.Opts.FixedRadius || size(frames, 1) < 3
            b = obj.Opts.Radius + b;
            ytemp = y(y>b & y<rows-b);
            xtemp = x(y>b & y<rows-b);
            y = ytemp(xtemp>b & xtemp<cols-b);
            x = xtemp(xtemp>b & xtemp<cols-b);
            r = ones(size(x)).*obj.Opts.Radius;
        else
            r = frames(3,:) .* obj.Opts.FrameScalingFactor;
            idx = zeros(size(x),'int8');
            for i=1:length(x)
                if r(i) < obj.Opts.Radius/2
                   r(i) = obj.Opts.Radius;
                end
                if y(i)-r(i)>b && y(i)+r(i)<rows-b && x(i)-r(i)>b && x(i)+r(i)<cols-b
                    idx(i) = 1;
                else
                    r(i) = min([y(i)-b-1, rows-b-y(i)-1, x(i)-b-1, cols-b-x(i)-1]);
                    if r(i) > obj.Opts.Radius/2
                        idx(i) = 1;
                    end
                end
            end
            y=y(idx>0);
            x=x(idx>0);
            r=double(r(idx>0));
        end
    end
  end
end