% LOCAL LEARNING OF 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 CsharedSketch.c;  % learning from aligned images 
mex CgetSUM2.c; % code for finding nearest neighbors, without inhibition, so very fast
mex CgetMAX1.c; % local maximum pooling for MAX1 maps
mex Csigmoid.c; % sigmoid transformation 
mex Ccopy.c; % copy around detect location
epsilon = .1; % allowed correlation between selected Gabors 
subsample = 1; % subsample in computing MAX1 and SUM2 maps   
locationShiftLimit = 2; % shift in normal direction = locationShiftLimit*subsample pixels
orientShiftLimit = 1; % shift in orientation
Correlation = CorrFilter(allFilter, epsilon); % correlation between filters 
numNeighbor = 5; % number of nearest neighbors including self
numElement = 60; % number of Gabors in active basis   
numIteration = 3;  % number of iterations
%% Load in training images
sizey = 120; % forced horizontal size of input image
imageFolder = 'positiveImage'; % folder of training images  
imageName = dir([imageFolder '\*.jpg']);
numImage = size(imageName, 1); % number of training images 
I = cell(1, numImage); 
maxSizex = 0; 
for (img = 1 : numImage)
    tmpIm = imread([imageFolder '\' imageName(img).name]); 
    if size(tmpIm,3) == 3
        tmpIm = rgb2gray(tmpIm);
    end
    [sx, sy] = size(tmpIm); 
    sizex = round(sx*sizey/sy); 
    if (maxSizex<sizex)
        maxSizex = sizex; 
    end
    I{img} = imresize(single(tmpIm), [sizex sizey], 'nearest'); 
end
sizex = maxSizex; 
allSizex = zeros(1, numImage) + sizex; 
allSizey = zeros(1, numImage) + sizey; 
disp(['number of training images = ' num2str(numImage)]);
%% Compute SUM1 maps and MAX1 maps
disp(['start filtering']); tic
[SUM1map, Ifill] = ApplyFilterfft(I, allFilter); % SUM1 maps by Gabor filtering
disp(['filtering time: ' num2str(toc) ' seconds']);
Csigmoid(numImage, allSizex, allSizey, numOrient, saturation, SUM1map);
MAX1map = cell(numImage, numOrient);
for (img = 1:numImage)
    for (orient = 1:numOrient) 
        MAX1map{img, orient} = single(zeros(floor(sizex/subsample), floor(sizey/subsample)));
    end
end
CgetMAX1(numImage, allSizex, allSizey, numOrient, locationShiftLimit, orientShiftLimit, subsample, ...
         SUM1map, MAX1map);           
SUM1mapLearn = cell(numNeighbor, numOrient); 
for (img = 1:numNeighbor)
    for (orient = 1:numOrient)
        SUM1mapLearn{img, orient} = single(zeros(sizex, sizey)); 
    end
end
%% 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, numNeighbor); % templates for training images 
for (img = 1 : numNeighbor)
    deformedTemplate{img} = single(zeros(sizex, sizey));  
end
SUM2score = zeros(numImage, 1); % template matching scores for training images 
allSUM2PositiveSelf = zeros(numImage, numImage); % matching score starting from each image
meanSUM2score = zeros(numImage, 1); % average matching score starting from each image
neighbor = zeros(numImage, numNeighbor); % neighbors of each starting image
allSelectedOrient = zeros(numImage, numElement);  % orientation and location of selected Gabors
allSelectedx = zeros(numImage, numElement); 
allSelectedy = zeros(numImage, numElement); 
allSelectedlambda = zeros(numImage, numElement); % weighting parameter for scoring template matching
allSelectedLogZ = zeros(numImage, numElement); % normalizing constant
%% Iteration part 0: learning template from this image
locationShiftLimit0 = 0; orientShiftLimit0 = 0; 
numImage0 = 1; 
SUM1map0 = cell(1, numOrient); 
for (orient = 1 : numOrient)
       SUM1map0{1, orient} = single(zeros(sizex, sizey));
end
showon = -1; 
for (this = 1:numImage)
ABlearnThis; 
end  % end of this 
%% Output locally learned templates with inhibition 
numInhibit = 5; OutTemplateInhibit; 
showon = 1; showNum = 0; 
for (this = selected)
    showNum = showNum + 1; 
    ABlearnThis; 
end
% selecting templates according to meanSUM2score. Once a template is selected, it inhibits the first numInhibit neighbors, 
% so that templates learned from these neighbors will not be selected. 
save 'locallyLearnedTemplate' numImage allTemplate allSUM2PositiveSelf ...
     allSelectedOrient allSelectedx allSelectedy allSelectedlambda allSelectedLogZ ...
     subsample locationShiftLimit orientShiftLimit numElement;  

