% 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
mex Csigmoid.c; % local normalization
mex CgetMAX1track.c; % local maximum pooling for MAX1 maps
epsilon = .1; % allowed correlation between selected Gabors 
subsample = 1; % subsample in computing MAX1 maps   
locationShiftLimit = 4; % shift in normal direction = locationShiftLimit*subsample pixels
orientShiftLimit = 1; % shift in orientation
numElement = 40; % number of Gabors in active basis   
Correlation = CorrFilter(allFilter, epsilon); % correlation between filters 
outputFolder = 'outputResult';
%% Load in training images
imageOption = 1; resizeFactor = .4; % resize images to specified ratio, and use the common upper-left portion 
%imageOption = 2; sizex = 100; sizey = 100; % resize images to specified common height and width
%imageOption = 3; sizey = 120; % resize images to have common width without changing their aspect ratios, and use the common horizontal central line
imageFolder = 'positiveImage'; % folder of training images  
imageName = dir([imageFolder '/*.jpg']);
numImage = size(imageName, 1); % number of training images 
I = cell(1, numImage); 
LoadImage; 
allSizex = zeros(1, numImage)+sizex; 
allSizey = zeros(1, numImage)+sizey; 
sizexSubsample = floor(sizex/subsample); 
sizeySubsample = floor(sizey/subsample); 
%% Compute SUM1 maps by Gabor filtering
disp(['start filtering']); tic
if (imageOption == 3) 
    [SUM1map, Ifill] = ApplyFilterfftFill(I, allFilter, localHalfx, localHalfy, thresholdFactor); % SUM1 maps by Gabor filtering
    I = Ifill; 
else 
    SUM1map = ApplyFilterfftSame(I, allFilter, localHalfx, localHalfy, -1, thresholdFactor);% SUM1 maps by Gabor filtering
end
Csigmoid(numImage, allSizex, allSizey, numOrient, saturation, SUM1map);
MAX1map = cell(numImage, numOrient);
trackMap = cell(numImage, numOrient);
for (img = 1:numImage)
    for (orient = 1:numOrient) 
        MAX1map{img, orient} = single(zeros(sizexSubsample, sizeySubsample));
        trackMap{img, orient} = single(zeros(sizexSubsample, sizeySubsample));
    end
end
CgetMAX1track(numImage, allSizex, allSizey, numOrient, locationShiftLimit, orientShiftLimit, subsample, ...
         SUM1map, MAX1map, trackMap); 
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, MAX1map, trackMap, ... % 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']);
%% Display and save results
OutLearningResult; 
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
%% Change aspect ratio of the template
allAspect = [.8 .9 1 1.2 1.4]; numAspect = 5; aspectOriginal = 3; 
allSelectedx = zeros(numAspect, numElement); 
allSelectedy = zeros(numAspect, numElement); 
allSelectedOrient = zeros(numAspect, numElement); 
for (a = 1:numAspect) 
   squeeze = allAspect(a); 
   if (a == aspectOriginal)
        for (i = 1 : numElement)
           allSelectedx(a, i) = selectedx(i); 
           allSelectedy(a, i) = selectedy(i); 
           allSelectedOrient(a, i) = selectedOrient(i); 
        end
   else 
        for (i = 1 : numElement)
           allSelectedx(a, i) = selectedx(i); 
           allSelectedy(a, i) = selectedy(i)*squeeze; 
           alpha = atan(tan(selectedOrient(i)*pi/numOrient)/squeeze); 
           if (alpha<0)
              alpha = alpha + pi; 
           end
           allSelectedOrient(a, i) = floor((alpha+.01)/pi*numOrient); 
        end
   end
end
save 'learnedTemplate' outputFolder numElement allSelectedOrient allSelectedx allSelectedy selectedlambda selectedLogZ ...
      numOrient halfFilterSize allFilter allSymbol locationShiftLimit orientShiftLimit saturation subsample ...
      localHalfx localHalfy thresholdFactor ...
      numAspect allAspect

