%% Compute SUM1 and MAX1 maps
close all; 
disp(['===============> Getting the SUM1 and MAX1 maps for testing images']);
for (img=1:numImage)
multipleResolutionImageName = ['working/multipleResolutionImage' num2str(img)];  
load(multipleResolutionImageName); 
for (s = 1 : numScale)
storeExponentialModelName = ['working/storedExponentialModel' num2str(s)];   
load(storeExponentialModelName); 
%% Prepare spaces for variables
allSizex = zeros(1, numResolution); allSizey = zeros(1, numResolution); 
MAX1maps = cell(numResolution, numOrient); 
trackMap = cell(numResolution, numOrient);
SUM2map = cell(1, numResolution); 
for(r=1:numResolution)
    [sizex, sizey] = size(J{r}); 
    allSizex(r) = sizex; allSizey(r) = sizey; % sizes of images at multiple resolutions
    sizexSubsample = floor(sizex/subsample); sizeySubsample = floor(sizey/subsample); 
    SUM2map{r} = single(zeros(floor((sizex+sizeTemplatex)/subsample), floor((sizey+sizeTemplatey)/subsample)));  % log-likelihood score maps for testing images
    averageMap{r} = single(zeros(sizex, sizey));
    for (orient=1:numOrient)
        MAX1map{r, orient} = single(zeros(sizexSubsample, sizeySubsample)); 
        trackMap{r, orient} = single(zeros(sizexSubsample, sizeySubsample));
    end
    translatedTemplate{r} = single(zeros(sizex, sizey));  % symbolic sketch on testing image
end
MAX2score = single(zeros(1, numResolution));  % maximum log-likelihood score at each resolution
Fx = zeros(1, numResolution); Fy = zeros(1, numResolution); % maximum likelihood location at each resolution
%% Compute SUM1 maps
disp(['start filtering image ' num2str(img) ' at Gabor length ' num2str(halfFilterSize*2+1) ' at all resolutions']); tic
SUM1map = ApplyFilterfft(J, allFilter, localOrNot, localHalfx, localHalfy, thresholdFactor);  % filter testing image at all resolutions
disp(['filtering time: ' num2str(toc) ' seconds']);
%% Compute average maps
Caverage(numResolution, allSizex, allSizey, numOrient, halfFilterSize, SUM1map, averageMap, halfTemplatex, halfTemplatey, windowNormalizeOrNot, thresholdFactor);
%% Compute MAX1 maps and track maps
disp(['start maxing image ' num2str(img) ' at Gabor length ' num2str(halfFilterSize*2+1) ' at all resolutions']); tic
CgetMAX1track(numResolution, allSizex, allSizey, numOrient, locationShiftLimit, orientShiftLimit, subsample, ...
              SUM1map, MAX1map, trackMap); 
disp(['maxing time: ' num2str(toc) ' seconds']);
%% Save the maps and spaces
SUM1MAX1mapName = ['working/SUM1MAX1map' 'image' num2str(img) 'scale' num2str(s)];   
save(SUM1MAX1mapName, 'SUM1map', 'MAX1map', 'trackMap', 'averageMap', 'SUM2map', 'translatedTemplate', 'MAX2score', 'Fx', 'Fy', 'allSizex', 'allSizey'); 
end
end
