%% Accumulate the SUM2 maps and compute MAX2 scores
clear all; close all; 
disp(['===============> Getting the SUM2 maps and MAX2 scores for testing images']);
load working/parameter; 
load working/testing; 
for (ii = 1:nn)
imageName = ['working/testingimage' num2str(ii)]; 
load(imageName); 
for (s = 1 : numScale)
  storeName = ['working/storedExponentialModel' num2str(s)];   
  load(storeName); 
  templateName = ['working/learnedTemplate' num2str(s)];   
  load(templateName); 
  mapName = ['working/SUMMAXmap' 'image' num2str(ii) 'scale' num2str(s)];  
  if (s==1) 
     load(mapName, 'SUM1map', 'MAX1map', 'trackMap', 'averageMap', 'SUM2map', 'MAX2score', 'Fx', 'Fy', 'allSizex', 'allSizey');
  else 
     load(mapName, 'SUM1map', 'MAX1map', 'trackMap', 'averageMap'); 
  end
disp(['start finding in image ' num2str(ii) ' at Gabor scale ' num2str(s) ' at all resolutions']); tic
CgetSUM2(numResolution, allSizex, allSizey, numOrient, subsample, ...
         sizeTemplatex, sizeTemplatey, numElement, selectedOrient, selectedx, selectedy, selectedlambda, selectedLogZ, ...
         MAX1map, SUM2map, averageMap, saturation, thresholdFactor);
disp(['finding time: ' num2str(toc) ' seconds']);
end
%% compute MAX2 scores for all resolutions
CgetMAX2(numResolution, allSizex, allSizey, subsample, ...
         sizeTemplatex, sizeTemplatey, SUM2map, MAX2score, Fx, Fy);
%% if inhibit, recalculate MAX2 score for each resolution
if (inhibitFind>0.)
    MAX2score = single(zeros(1, numResolution));  % recalculate MAX2 score at each resolution
   for (s = 1 : numScale)
      storeName = ['working/storedExponentialModel' num2str(s)];   
      load(storeName); 
      templateName = ['working/learnedTemplate' num2str(s)];   
      load(templateName); 
      mapName = ['working/SUMMAXmap' 'image' num2str(ii) 'scale' num2str(s)];      
      load(mapName, 'SUM1map', 'MAX1map', 'trackMap', 'averageMap'); 
      CgetBackInhibit(numResolution, allSizex, allSizey, numOrient, ...
         locationShiftLimit, orientShiftLimit, subsample, halfFilterSize, allSymbol(1, :), ...
         numElement, selectedOrient, selectedx, selectedy, selectedlambda, selectedLogZ, ...
         SUM1map, MAX1map, trackMap, translatedTemplate, Fx, Fy, saturation, ...
         sizeTemplatex, sizeTemplatey, Correlation, MAX2score, averageMap);
   end
end
%% find the optimal resolution
[maxOverResolution ind] = max(MAX2score);   % most likely resolution
Outepsgif(J{ind},  ['output/testImage' num2str(ii) 'original' num2str(s)]); % original image 
%% plot SUM2 map and MAX2 scores
figure('Position', [1 1 1024 768]); 
surf(double(SUM2map{ind})); % log-likelihood score map at optimal resolution
set(gca,'FontSize',18); % adjust the display format and save the plot
set(gcf, 'PaperPositionMode', 'auto');
print('-depsc', ['output/testImage' num2str(ii) 'SUM2mapUsingGaborScale' num2str(s) '.eps']);
figure; h1=plot(MAX2score); 
title('MAX2 score vs scale'); xlabel('scale'); ylabel('MAX2');
set(gca,'FontSize',16); % adjust the display format and save the plot
set(get(gca,'XLabel'),'FontSize',18);
set(get(gca,'YLabel'),'FontSize',18);
set(get(gca,'Title'),'FontSize',20);
set(h1,'Linewidth',3);
set(gcf, 'PaperPositionMode', 'auto');
print('-depsc', ['output/testImage' num2str(ii) 'MAX2scoreUsingGaborScale' num2str(s) '.eps']);
%% superpose templates at multiple scales
if (inhibitFind>0.)
    MAX2score = single(zeros(1, numResolution));  % recalculate MAX2 score at each resolution
end
for (s = 1 : numScale)
      storeName = ['working/storedExponentialModel' num2str(s)];   
      load(storeName); 
      templateName = ['working/learnedTemplate' num2str(s)];   
      load(templateName); 
      mapName = ['working/SUMMAXmap' 'image' num2str(ii) 'scale' num2str(s)];       
      load(mapName, 'SUM1map', 'MAX1map', 'trackMap', 'averageMap', 'translatedTemplate'); 
      if (inhibitFind>0.)
         CgetBackInhibit(numResolution, allSizex, allSizey, numOrient, ...
         locationShiftLimit, orientShiftLimit, subsample, halfFilterSize, allSymbol(1, :), ...
         numElement, selectedOrient, selectedx, selectedy, selectedlambda, selectedLogZ, ...
         SUM1map, MAX1map, trackMap, translatedTemplate, Fx, Fy, saturation, ...
         sizeTemplatex, sizeTemplatey, Correlation, MAX2score, averageMap);
      else
         CgetBack(numResolution, allSizex, allSizey, numOrient, ...
         locationShiftLimit, orientShiftLimit, subsample, halfFilterSize, allSymbol(1, :), ...
         numElement, selectedOrient, selectedx, selectedy, ...
         MAX1map, trackMap, translatedTemplate, Fx, Fy, saturation);
      end
      Outepsgif(J{ind}+translatedTemplate{ind}*100, ['output/testImage' num2str(ii) 'DetectionUsingGaborScale' num2str(s)]); % overlaid sketch at optimal resolution  
end
end

