%% Accumulate the SUM2 maps and compute MAX2 scores
disp(['===============> Getting the SUM2 maps and MAX2 scores for testing images']);
for (img = 1:numImage)
close all;
if (img~=starting)    
multipleResolutionImageName = ['working/multipleResolutionImage' num2str(img)];
load(multipleResolutionImageName); 
for (s = 1 : numScale)
  storeExponentialModelName = ['working/storedExponentialModel' num2str(s)];   
  load(storeExponentialModelName); 
  learnedTemplateName = ['working/learnedTemplate' num2str(s)];   
  load(learnedTemplateName); 
  SUM1MAX1mapName = ['working/SUM1MAX1map' 'image' num2str(img) 'scale' num2str(s)]; 
  if (s==1) 
     load(SUM1MAX1mapName, 'SUM1map', 'MAX1map', 'trackMap', 'averageMap', 'SUM2map', 'MAX2score', 'Fx', 'Fy', 'allSizex', 'allSizey');
  else 
     load(SUM1MAX1mapName, 'SUM1map', 'MAX1map', 'trackMap', 'averageMap'); 
  end
disp(['start finding in image ' num2str(img) ' at Gabor length ' num2str(halfFilterSize*2+1) ' at all resolutions']); tic
CgetSUM2(numResolution, allSizex, allSizey, numOrient, subsample, ...
         sizeTemplatex, sizeTemplatey, numElement, selectedOrient, selectedx, selectedy, selectedlambda, selectedLogZ, ...
         MAX1map, SUM2map, averageMap, saturation);
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)
      storeExponentialModelName = ['working/storedExponentialModel' num2str(s)];   
      load(storeExponentialModelName); 
      learnedTemplateName = ['working/learnedTemplate' num2str(s)];   
      load(learnedTemplateName); 
      SUM1MAX1mapName = ['working/SUM1MAX1map' 'image' num2str(img) 'scale' num2str(s)];       
      load(SUM1MAX1mapName, '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
%% plot SUM2 map and MAX2 scores
if (showOriginal>0.)
Outepsgif(J{ind},  ['output/Image' num2str(img) 'original' num2str(s)]); % original image 
end
if (showMAX2SUM2>0.)
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/Image' num2str(img) '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/Image' num2str(img) 'MAX2scoreUsingGaborScale' num2str(s) '.eps']);
end
%% superpose templates at multiple scales
if (inhibitFind>0.)
    MAX2score = single(zeros(1, numResolution));  % recalculate MAX2 score at each resolution
end
for (s = 1 : numScale)
      storeExponentialModelName = ['working/storedExponentialModel' num2str(s)];   
      load(storeExponentialModelName); 
      learnedTemplateName = ['working/learnedTemplate' num2str(s)];   
      load(learnedTemplateName); 
      SUM1MAX1mapName = ['working/SUM1MAX1map' 'image' num2str(img) 'scale' num2str(s)];     
      load(SUM1MAX1mapName, '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/LNlion' num2str(img)]); % overlaid sketch at optimal resolution  
%% Crop detected window for learning
learningSpaceName = ['working/learningSpace' num2str(s)]; 
load(learningSpaceName); 
maxFx = Fx(ind); maxFy = Fy(ind); 
sizex = allSizex(ind); sizey = allSizey(ind); 
ave = averageMap{ind}; 
xHere = maxFx*subsample+halfTemplatex+1; yHere = maxFy*subsample+halfTemplatey+1; 
if ((xHere>=1)&&(xHere<=sizex)&&(yHere>=1)&&(yHere<=sizey))
     averageHere = ave(xHere, yHere);  % average within local window
else
     averageHere = single(1.); 
end
for (orient = 1 : numOrient)
   Ccopy(SUM1mapLearn{img, orient}, SUM1map{ind, orient}, maxFx*subsample, maxFy*subsample, sizeTemplatex, sizeTemplatey, sizex, sizey,  averageHere, saturation); 
  % Ccopy(MAX1mapLearn{img, orient}, MAX1map{ind, orient}, maxFx, maxFy, sizeTemplatexSubsample, sizeTemplateySubsample, floor(sizex/subsample), floor(sizey/subsample),  averageHere, saturation); 
  % Ccopy0(trackMapLearn{img, orient}, trackMap{ind, orient}, maxFx, maxFy, sizeTemplatexSubsample, sizeTemplateySubsample, floor(sizex/subsample), floor(sizey/subsample)); 
end
save(learningSpaceName, 'SUM1mapLearn', 'MAX1mapLearn', 'trackMapLearn'); 
end
end
end


