% SHARED SKETCH ALGORITHM FOR LEARNING ACTIVE BASIS
%% Load in exponential model, mex C codes, and set parameters
% ExponentialModel; % this line is not needed if 'storedExponentialModel' exists
clear; close all; 
load 'storedExponentialModel'; % load in exponential model 
mex Clearn.c; % learning by shared sketch algorithm
epsilon = .1; % allowed correlation between selected Gabors 
subsample = 2; % subsample in computing MAX1 maps   
locationShiftLimit = 3; % shift in normal direction = locationShiftLimit*subsample pixels
orientShiftLimit = 1; % shift in orientation
numElement = 50; % number of Gabors in active basis   
Correlation = CorrFilter(allFilter, epsilon); % correlation between filters 
outputFolder = 'butterflyEXP3';
mkdir([outputFolder 'eps']); 
mkdir([outputFolder 'png']); 
%% Load in training images
useMinimumSize = 1; resizeFactor = 1.; % resize the input image 
useFixedSize = -1; sizex = 82; sizey = 164; % forced size of inpu image
imageFolder = 'positiveImage'; % folder of training images  
imageName = dir([imageFolder '/*.jpg']);
numImage = size(imageName, 1); % number of training images 
I = cell(1, numImage); 
for (img = 1 : numImage)
    tmpIm = imread([imageFolder '/' imageName(img).name]); 
    if size(tmpIm,3) == 3
        tmpIm = rgb2gray(tmpIm);
    end
    if (useMinimumSize>0)
       I{img} = imresize(single(tmpIm), resizeFactor, 'nearest');
    end
    if (useFixedSize>0)
       I{img} = imresize(single(tmpIm), [sizex sizey], 'nearest');      
    end     
end
if (useMinimumSize>0)   
    imageSize = zeros(numImage, 2);
    for (img = 1 : numImage)
      imageSize(img, :) = size(I{img});   
    end
    sizex = min(imageSize(:, 1)); sizey = min(imageSize(:, 2)); 
    for (img = 1 : numImage)
      I{img} = I{img}(1:sizex, 1:sizey); 
    end
end  

%% Compute SUM1 maps by Gabor filtering
disp(['start filtering']); tic
SUM1map = ApplyFilterfftSame(I, allFilter, localOrNot, localHalfx, localHalfy, -1, thresholdFactor);% SUM1 maps by Gabor filtering
disp(['filtering time: ' num2str(toc) ' seconds']);
%% Prepare output variables for learning
selectedOrient = zeros(1, numElement);  % orientation and location of selected Gabors
selectedx = zeros(1, numElement); 
selectedy = zeros(1, numElement); 
selectedlambda = zeros(1, numElement); % weighting parameter for scoring template matching
selectedLogZ = zeros(1, numElement); % normalizing constant
commonTemplate = single(zeros(sizex, sizey)); % template of active basis 
deformedTemplate = cell(1, numImage); % templates for training images 
for (img = 1 : numImage)
    deformedTemplate{img} = single(zeros(sizex, sizey));  
end
SUM2score = zeros(numImage, 1); % template matching scores for training images 
%% Learning by shared sketch algorithm
disp(['start mex-C learning']); tic
Clearn(numOrient, locationShiftLimit, orientShiftLimit, saturation, subsample, ... % about active basis  
       numElement, numImage, sizex, sizey, SUM1map, ... % about training images 
       halfFilterSize, Correlation, allSymbol(1, :), ... % about filters
       numStoredPoint, storedlambda, storedExpectation, storedLogZ, ... % about exponential model 
       selectedOrient, selectedx, selectedy, selectedlambda, selectedLogZ, SUM2score, ... % learned parameters
       commonTemplate, deformedTemplate); % learned templates 
disp(['mex-C learning time: ' num2str(toc) ' seconds']);
%% save results in .m
minx = min(selectedx); maxx = max(selectedx); centerx = floor((minx+maxx)/2.); 
miny = min(selectedy); maxy = max(selectedy); centery = floor((miny+maxy)/2.); 
for (i = 1 : numElement)
    selectedx(i) = selectedx(i) - centerx; 
    selectedy(i) = selectedy(i) - centery; 
end
save 'learnedTemplate' numElement selectedOrient selectedx selectedy selectedlambda selectedLogZ ...
      numOrient halfFilterSize allFilter allSymbol locationShiftLimit orientShiftLimit saturation subsample ...
      localOrNot localHalfx localHalfy thresholdFactor outputFolder; 
%% save figures into eps and png folders  
pic = -double(commonTemplate); showImage; 
saveas(gcf, [outputFolder 'eps/template.eps'], 'eps');
saveas(gcf, [outputFolder 'png/template.png'], 'png');
pic = commonTemplate*0.+255; showImage; 
saveas(gcf, [outputFolder 'eps/empty.eps'], 'eps');
saveas(gcf, [outputFolder 'png/empty.png'], 'png');
for (img = 1 : numImage)
    close all; 
    tmpIm = imresize(imread([imageFolder '/' imageName(img).name]), resizeFactor, 'nearest');
    pic = tmpIm(1:sizex, 1:sizey, :); showImage; 
    saveas(gcf, [outputFolder 'eps/input' num2str(100+img) '.eps'], 'epsc');
    saveas(gcf, [outputFolder 'png/input' num2str(100+img) '.png'], 'png');    
    pic = -(deformedTemplate{img}); showImage; 
    saveas(gcf, [outputFolder 'eps/sketch' num2str(100+img) '.eps'], 'eps');
    saveas(gcf, [outputFolder 'png/sketch' num2str(100+img) '.png'], 'png');
end
%% write the html code for reproducibility page
heightstr = '"height=90> '; 
disp(['<IMG SRC="' outputFolder 'png/template' '.png' heightstr]);
disp(['<IMG SRC="' outputFolder 'png/empty' '.png' heightstr]);
counter = 1; length = 5; 
for (img = 1 : numImage)
 disp(['<IMG SRC="' outputFolder 'png/input' num2str(100+img) '.png' heightstr]);
 disp(['<IMG SRC="' outputFolder 'png/sketch' num2str(100+img) '.png' heightstr]);
 counter = counter + 1; 
 if (counter == length)
     disp('<BR>'); 
     counter = 0; 
 end
end
%% write the LaTex code for paper 
heightstr = 'height=.6in}'; 
disp(['\psfig{figure=' outputFolder 'eps/template' '.eps,'  heightstr]); 
disp(['\psfig{figure=' outputFolder 'eps/empty' '.eps,'  heightstr]);
counter = 1; length = 5; 
for (img = 1 : numImage)
 disp(['\psfig{figure=' outputFolder 'eps/input' num2str(100+img) '.eps,' heightstr]);
 disp(['\psfig{figure=' outputFolder 'eps/sketch' num2str(100+img) '.eps,' heightstr]);
 counter = counter + 1; 
 if (counter == length)
     disp('\\'); 
     counter = 0; 
 end
end


