%% 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,'J0');

    storeExponentialModelName = ['storedExponentialModel'];   
    load(storeExponentialModelName);
    %% Prepare spaces for variables
    
    %% Compute SUM1 maps
    % use non-inscribed images to compute SUM1 maps, to avoid boundary effects 
    disp(['start filtering image ' num2str(img) ' at Gabor length ' num2str(halfFilterSize*2+1) ' at all resolutions']); tic
    SUM1map = ApplyFilterfft({J0}, allFilter, localOrNot, localHalfx, localHalfy, thresholdFactor);  % filter testing image at all resolutions
    disp(['filtering time: ' num2str(toc) ' seconds']);
    %% inscribe the SUM1 maps into a square bounding box
    % change allSizex and allSizey

    for j = 1:size(SUM1map,2)
        currentOuterBBsize = max(size(SUM1map{1,j}));
        SUM1map{1,j} = single(inscribe([currentOuterBBsize currentOuterBBsize], SUM1map{1,j}, 0));
    end
   
    
    %% Still, prepare spaces for variables
    [sizex, sizey] = size(SUM1map{1});
    sizexSubsample = floor(sizex/subsample); sizeySubsample = floor(sizey/subsample);
    averageMap = single(zeros(sizex, sizey));
    for (orient=1:numOrient)
        MAX1map{1, orient} = single(zeros(sizexSubsample, sizeySubsample)); 
        trackMap{1, orient} = single(zeros(sizexSubsample, sizeySubsample));
    end
    translatedTemplate{1} = single(zeros(sizex, sizey));  % symbolic sketch on testing image
    
    %% Compute MAX1 maps and track maps
    disp(['start maxing image ' num2str(img) ' at Gabor length ' num2str(halfFilterSize*2+1) ' at all resolutions']); tic
    subsampleM1 = 1;  % to be safe, please refrain from setting it to be larger than 1 
	[MAX1map ARGMAX1map M1RowShift M1ColShift M1OriShifted] = ...
		mexc_ComputeMAX1( SUM1map, locationShiftLimit,...
        orientShiftLimit, subsampleM1 );
	
	% sigmoid transform on S1 and M1 maps
	mexc_Sigmoid( saturation, SUM1map );
	mexc_Sigmoid( saturation, MAX1map );

    disp(['maxing time: ' num2str(toc) ' seconds']);
    %% Save the maps and spaces
    SUM1MAX1mapName = ['working/SUM1MAX1map' 'image' num2str(img)];   
    save(SUM1MAX1mapName, 'SUM1map', 'MAX1map', 'M1RowShift', 'M1ColShift',...
		'M1OriShifted', 'ARGMAX1map', 'translatedTemplate');
end
