% SUM-MAX MAPS FOR TEMPLATE MATCHING INFERENCE
%% load in learned template, set parameter values
clear all; 
load 'learnedTemplate'; % load in the learned template
mex CgetSUM2.c;   % compile C code for template matching inference
mex CgetMAX1.c; 
%% load in test image, generate multiple resolutions
numResolution = 10;  % total number of resolutions to search 
allSizex = zeros(1, numResolution); allSizey = zeros(1, numResolution); 
imageFolder = 'testImage'; % folder of training images  
imageName = dir([imageFolder '/*.jpg']);
numImage = size(imageName, 1); % number of training images 
for (ii = 1 : numImage)
    close all; 
    tmpIm = imread([imageFolder '/' imageName(ii).name]); 
    pic = tmpIm; showImage;
saveas(gcf, [outputFolder 'eps/testoriginal' num2str(100+ii) '.eps'], 'epsc');
saveas(gcf, [outputFolder 'png/testoriginal' num2str(100+ii) '.png'], 'png');
    if size(tmpIm,3) == 3
        tmpIm = rgb2gray(tmpIm);
    end
J0 = single(tmpIm); 
MAX1maps = cell(numResolution, numOrient); 
SUM2map = cell(1, numResolution); 
for(img=1:numResolution)
    J{img} = imresize(J0, .10+(img-1)*.1, 'nearest');  % scaled images
    [sizex, sizey] = size(J{img}); 
    allSizex(img) = sizex; allSizey(img) = sizey; % sizes of images at multiple resolutions
    SUM2map{img} = single(zeros(sizex, sizey));  % log-likelihood score maps for testing images
    for (orient=1:numOrient)
        MAX1map{img, orient} = single(zeros(sizex, sizey)); 
    end
    translatedTemplate{img} = single(zeros(sizex, sizey));  % symbolic sketch on testing image
end
%% compute SUM1 maps for all resolutions
disp(['start filtering']); tic
SUM1map = ApplyFilterfft(J, allFilter, localHalfx, localHalfy, thresholdFactor); % filtering images at multiple resolutions
disp(['filtering time: ' num2str(toc) ' seconds']);
MAX2score = single(zeros(1, numResolution));  % maximum log-likelihood score at each resolution
%% detect the object and superpose template
disp(['start sigmoid and maxing']); tic
Csigmoid(numResolution, allSizex, allSizey, numOrient, saturation, SUM1map);
CgetMAX1(numResolution, allSizex, allSizey, numOrient,  ...
                     locationShiftLimit, orientShiftLimit, ...
                     SUM1map, MAX1map);
disp(['sigmoid and maxing time: ' num2str(toc) ' seconds']);     
disp(['start finding']); tic
MMAX2 = -1e10; 
for (flip = 0 : 1)
   for (rot = -rotateShiftLimit : rotateShiftLimit)
         r = rot+rotateShiftLimit+1 + (rotateShiftLimit*2+1)*flip;  
         MAX2score = single(zeros(1, numResolution)); Fx = zeros(1, numResolution); Fy = zeros(1, numResolution);      
         CgetSUM2(1, 1, numResolution, allSizex, allSizey, numOrient, ...
              locationShiftLimit, orientShiftLimit, ...
              sizeTemplatex, sizeTemplatey, ...
              halfFilterSize, SUM1map, allSymbol(1, :), ...
              numElement, allSelectedOrient(r, :), allSelectedx(r, :), allSelectedy(r, :), selectedlambda, selectedLogZ, ...
              MAX1map, SUM2map, translatedTemplate, MAX2score, Fx, Fy);
         [maxOverResolution j] = max(MAX2score);   % most likely resolution
         if (MMAX2 < maxOverResolution) 
            MMAX2 = maxOverResolution; Mrot = rot; Mflip = flip; 
         end          
disp(['mex-C finding time: ' num2str(toc) ' seconds']);
%% display results
         [maxOverResolution ind] = max(MAX2score);   % most likely resolution
         if (MMAX2 < maxOverResolution) 
               MMAX2 = maxOverResolution; Mrot = rot; 
         end
   end
end
        r = Mrot+rotateShiftLimit+1 + (rotateShiftLimit*2+1)*Mflip;  tic; 
        MAX2score = single(zeros(1, numResolution)); Fx = zeros(1, numResolution); Fy = zeros(1, numResolution); 
        CgetSUM2(1, 1, numResolution, allSizex, allSizey, numOrient, ...
              locationShiftLimit, orientShiftLimit, ...
              sizeTemplatex, sizeTemplatey, ...
              halfFilterSize, SUM1map, allSymbol(1, :), ...
              numElement, allSelectedOrient(r, :), allSelectedx(r, :), allSelectedy(r, :), selectedlambda, selectedLogZ, ...
              MAX1map, SUM2map, translatedTemplate, MAX2score, Fx, Fy);
        disp(['mex-C finding time: ' num2str(toc) ' seconds']);
[maxOverResolution ind] = max(MAX2score);   % most likely resolution
pic = double(J{ind}+translatedTemplate{ind}*100); showImage;
saveas(gcf, [outputFolder 'eps/testsketch' num2str(100+ii) '.eps'], 'eps');
saveas(gcf, [outputFolder 'png/testsketch' num2str(100+ii) '.png'], 'png');
if (0)
figure('Position', [1 1 1024 768]); 
mesh(double(SUM2map{ind})); % log-likelihood score map at optimal resolution
saveas(gcf, 'ex2_1_S2','eps');
figure; plot(MAX2score); 
title('MAX2 score vs resolution'); xlabel('resolution'); ylabel('MAX2');
saveas(gcf, 'ex2_1_M2','eps');
end
end

fid = fopen([outputFolder '.html'], 'at'); 
   outMessage =  ['==> Number of testing images = ' num2str(numImage) ...
                          '; Range of rotation = ' num2str(rotateShiftLimit) ' times pi/15'  ...
                          '<br>'];
    fprintf(fid, '%s\n', outMessage);
    disp(outMessage);           
%% write the html code for reproducibility page

heightstr = '"height=160> '; 
counter = 0; 
for (ii = 1 : numImage)
 fprintf(fid, '%s\n', ['<IMG SRC="' outputFolder 'png/testoriginal' num2str(100+ii) '.png' heightstr]);
 fprintf(fid, '%s\n', ['<IMG SRC="' outputFolder 'png/testsketch' num2str(100+ii) '.png' heightstr]);
 counter = counter + 1; 
 if (counter == 3)
     fprintf(fid, '%s\n', ['<br>']);
     counter = 0; 
 end
end
    fprintf(fid, '%s\n', ['<hr>']);
    fclose(fid);



