% clear all;
% close all;
% mex MatchingPursuit.c;

function [template, numSelectedFeatures] = OrthogonalFiltersSelection(I, GarborScaleList, DoGScaleList, numOrient, numSketch, locationShiftLimit, orientShiftLimit, sx, sy,isOrthogonal,subsamplingOMP, LocalNormOrNot, scaleFactor,thresholdFactor,FRAMEfilterBankInfo)

%% set parameters 

numScale = length(GarborScaleList); %number of scales for Gabor wavelets
numDoG = length(DoGScaleList);

epsilon = .1; 
Upperbound = 6.; 
Shutup = 0; 

numImg = length(I);

%% read training images
% imgFold = 'Image';
% imgName = dir([imgFold '/*.jpg']);
% numImg = size(imgName, 1);
% clear I;
% for (i = 1 : numImg)
%     img = imread([imgFold '/' imgName(i).name]);
%     img = imresize(img, [sx,sy],'nearest');
%     img = im2double(rgb2gray(img));
%     img=img-mean(img(:));
%     img=img/std(img(:));
%     I{i} = img;
% end


%% filter images
tic;
[allFilterR, allFilterI, allSymbol] = MakeFilterBank(GarborScaleList, DoGScaleList, numOrient);
disp(['finished generating filter bank: '  num2str(toc) ' seconds']);

tic;
[Crr, Cri, Cir, Cii] = CorrFilter(allFilterR, allFilterI); %calculate correlation among filters
disp(['finished computing filters correlation: ' num2str(toc) ' seconds']);

tic
CorrAll = CorrFilterFrame([allFilterR, allFilterI]); %calculate correlation among filters (This is for the orthogonal matching pursuit)
disp(['finished computing another filters correlation: ' num2str(toc) ' seconds']);

tic;
[allFilteredR, allFilteredI, allFiltered] = ApplyFilter(I, allFilterR, allFilterI);
disp(['finished images filtering: ' num2str(toc) ' seconds']);

half = zeros(1, numOrient*numScale+numDoG); 
for (i=1:numOrient*numScale+numDoG)
    half(i) = (size(allFilterR{1, i}, 1)-1)/2;
end



% local normalization
if LocalNormOrNot
    % 
    
    for iI = 1 : numImg
      for iScale=1:numScale  
        h = half(1+(iScale-1)*numOrient); 
        
        [allFilteredR(iI, 1+(iScale-1)*numOrient : numOrient+(iScale-1)*numOrient)]...
            = LocalNormalizeDouble(allFilteredR(iI, 1+(iScale-1)*numOrient : numOrient+(iScale-1)*numOrient),[],h,round(scaleFactor*h),round(scaleFactor*h),thresholdFactor);    
        
        [allFilteredI(iI, 1+(iScale-1)*numOrient : numOrient+(iScale-1)*numOrient)]...
            = LocalNormalizeDouble(allFilteredI(iI, 1+(iScale-1)*numOrient : numOrient+(iScale-1)*numOrient),[],h,round(scaleFactor*h),round(scaleFactor*h),thresholdFactor);
      end
    
      % for DoG
      for iScale=1:numDoG  
        ind  = numOrient*numScale+iScale;
        h = half(ind);       
        [allFilteredR(iI, ind)]=LocalNormalizeDouble(allFilteredR(iI, ind),[],h,round(scaleFactor*h),round(scaleFactor*h),thresholdFactor);                
      end
    end
    
    for iI = 1 : numImg
        for iF = 1 : numOrient*numScale+numDoG
             allFiltered{iI,iF}=(allFilteredR{iI,iF}.^2 + allFilteredI{iI,iF}.^2)/2;
        end
    end
 
end


%% prepare for storing experimental results
% multiple-scaled Gobar template
syma = zeros(sx, sy); 

if numDoG>0
    numDoGTemplate=1;
else
    numDoGTemplate=0;
end

sym = cell(1, numScale+numDoGTemplate); 
for (i = 1 : numScale+numDoGTemplate)
    sym{i} = zeros(sx, sy); 
end

% recontructed images
Asym = cell(1, numImg); 
for (i = 1 : numImg)
    Asym{i} = zeros(sx, sy); 
end

% deformed Gabor templates
Asyma = cell(1, numImg); 
for (i = 1 : numImg)
    Asyma{i} = zeros(sx, sy); 
end

% deformed DoG templates
Asym0 = cell(1, numImg*(numScale+1));
for (i = 1 : numImg*(numScale+1))
    Asym0{i} = zeros(sx, sy);  
end

filterSelected = zeros(1, numSketch);
xSelected = zeros(1, numSketch);
ySelected = zeros(1, numSketch);


%% do matching pursuit
tic;
% MatchingPursuit(numImg, numOrient*numScale+numDoG, numOrient, numDoG, allFilteredR,...
%                 allFilteredI, allFiltered, Crr, Cri, Cir, Cii, half, sx, sy,...
%                 allFilterR(1, :), allFilterI(1, :), allSymbol(1, :), syma, sym,...
%                 Asym, Asym0, Asyma, locationShiftLimit, orientShiftLimit, numSketch,...
%                 epsilon, Upperbound, SHUTUP, filterSelected, xSelected, ySelected, numSketchSelected);


deformed_filterSelected=zeros(numImg,numSketch);
deformed_xSelected=zeros(numImg,numSketch);
deformed_ySelected=zeros(numImg,numSketch);
coefficientsSinList=zeros(numImg,numSketch);
coefficientsCosList=zeros(numImg,numSketch);

% if (isOrthogonal == false)
%     disp(['Normal Matching Pursuit']);
%     
% numSketchSelected = MatchingPursuit(numImg, numOrient*numScale+numDoG, numOrient, numDoG, allFilteredR,...
%                 allFilteredI, allFiltered, Crr, Cri, Cir, Cii, half, sx, sy,...
%                 allFilterR(1, :), allFilterI(1, :), allSymbol(1, :), syma, sym,...
%                 Asym, Asym0, Asyma, locationShiftLimit, orientShiftLimit, numSketch,...
%                 epsilon, Upperbound, Shutup, filterSelected, xSelected, ySelected);
%            
% else
%     disp(['Orthogonal Matching Pursuit']);
%    
% numSketchSelected = OrthogonalMatchingPursuit(I, numOrient*numScale+numDoG, numOrient, numDoG, allFilteredR,...
%                 allFilteredI, allFiltered, Crr, Cri, Cir, Cii, half, sx, sy,...
%                 allFilterR(1, :), allFilterI(1, :), allSymbol(1, :), syma, sym,...
%                 Asym, Asym0, Asyma, locationShiftLimit, orientShiftLimit, numSketch,...
%                 epsilon, Shutup,Upperbound,  filterSelected, xSelected, ySelected, ...
%                 deformed_filterSelected, deformed_xSelected, deformed_ySelected, ...
%                 coefficientsSinList, coefficientsCosList, CorrAll, subsamplingOMP);      
%    
% end 

if (isOrthogonal == false)
    disp(['Normal Matching Pursuit']);
    actualSubsamplingOMP = numSketch+1;     % when the SubsamplingOMP is larger than the number of filters we want to selelct, then OMP performs like MP.
else
    disp(['Orthogonal Matching Pursuit with subsample ' num2str(subsamplingOMP)]);
    actualSubsamplingOMP = subsamplingOMP;
end
  
numSketchSelected = OrthogonalMatchingPursuit(I, numOrient*numScale+numDoG, numOrient, numDoG, allFilteredR,...
                allFilteredI, allFiltered, Crr, Cri, Cir, Cii, half, sx, sy,...
                allFilterR(1, :), allFilterI(1, :), allSymbol(1, :), syma, sym,...
                Asym, Asym0, Asyma, locationShiftLimit, orientShiftLimit, numSketch,...
                epsilon, Shutup,Upperbound,  filterSelected, xSelected, ySelected, ...
                deformed_filterSelected, deformed_xSelected, deformed_ySelected, ...
                coefficientsSinList, coefficientsCosList, CorrAll, actualSubsamplingOMP);      

disp(['finished Matching Pursuit: ' num2str(toc) ' seconds']);
  
            
 filterSelected = filterSelected(1:numSketchSelected); 
 xSelected = xSelected(1:numSketchSelected);
 ySelected = ySelected(1:numSketchSelected);

 %[indEliminated, numSelectedFeatures]= createMask(filterSelected,xSelected,ySelected, sx, sy, numOrient, numScale, numDoG);
 [filterSelected, xSelected,ySelected, numSelectedFeatures]= indexTtransformation(filterSelected, xSelected, ySelected, sx, sy, numOrient, numScale, numDoG);

 template=struct('selectedFilter',single(filterSelected),'selectedRow',single(xSelected),'selectedCol',single(ySelected));
 
 
 if LocalNormOrNot
     
    ReconOriginal=cell(1,numImg);
    
    Corr = FRAMEfilterBankInfo.Corr;
    halfFilterSize = FRAMEfilterBankInfo.halfFilterSize;
    filters = FRAMEfilterBankInfo.filters;
    
    deformed_filterSelected = deformed_filterSelected(:,1:numSketchSelected); 
    deformed_xSelected = deformed_xSelected(:,1:numSketchSelected);
    deformed_ySelected = deformed_ySelected(:,1:numSketchSelected);
    
    for (iImg=1:numImg)
      [deformedfilter, deformedx,deformedy, ~]= indexTtransformation(deformed_filterSelected(iImg,:), ...
          deformed_xSelected(iImg,:), deformed_ySelected(iImg,:), sx, sy, numOrient, numScale, numDoG);
      temp=struct('selectedFilter',single(deformedfilter),'selectedRow',single(deformedx),'selectedCol',single(deformedy));
      [~, ~, ~, ReconOriginal{iImg}]=ComputeProjection(Corr,temp,single(halfFilterSize),filters,single(I{iImg}));
    end
    
     
 end
 
 
%% show and save experimental results
ShowAndSaveResultsV2;

%% generate html to show filter selection experimental results
GenerateHtmlV2;




