% 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 Clearn.c;  % learning from aligned images 
mex Cmatch.c; % code for finding nearest neighbors, without inhibition, so very fast
mex CsigmoidTransform.c; % sigmoid transformation 
mex ClocalMaximumPooling.c; % local maximum pooling for MAX1 maps
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 = 40; % number of Gabors in active basis   
numIteration = 3;  % number of iterations
%% Load in training images
useMinimumSize = -1; resizeFactor = 1.; % resize the input image 
useFixedSize = 1; sizex = 100; sizey = 100; % 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        
allSizex = zeros(1, numImage) + sizex; 
allSizey = zeros(1, numImage) + sizey; 
sizexSubsample = floor(sizex/subsample); 
sizeySubsample = floor(sizey/subsample); 
disp(['number of training images = ' num2str(numImage)]);
%% Compute SUM1 maps and MAX1 maps
disp(['start filtering']); tic
SUM1map = ApplyFilterfftSame(I, allFilter); % SUM1 maps by Gabor filtering
disp(['filtering time: ' num2str(toc) ' seconds']);
CsigmoidTransform(numImage, allSizex, allSizey, numOrient, saturation, SUM1map);
MAX1map = cell(numImage, numOrient);
tracking = cell(numImage, numOrient); 
for (img = 1:numImage)
    for (orient = 1:numOrient) 
        MAX1map{img, orient} = single(zeros(sizexSubsample, sizeySubsample));
        tracking{img, orient} = single(zeros(sizexSubsample, sizeySubsample));
    end
end
ClocalMaximumPooling(numImage, allSizex, allSizey, numOrient,  ...
                     locationShiftLimit, orientShiftLimit, subsample, saturation, ...
                     SUM1map, MAX1map, tracking);           
SUM1mapLearn = cell(numNeighbor, numOrient); 
MAX1mapLearn = cell(numNeighbor, numOrient); 
trackingLearn = cell(numNeighbor, numOrient); 
for (img = 1:numNeighbor)
    for (orient = 1:numOrient)
        SUM1mapLearn{img, orient} = single(zeros(sizex, sizey)); 
        MAX1mapLearn{img, orient} = single(zeros(sizexSubsample, sizeySubsample)); 
        trackingLearn{img, orient} = single(zeros(sizexSubsample, sizeySubsample)); 
    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 
allSUM2score = 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
%% Iteration part 0: learning template from this image
locationShiftLimit0 = 0; orientShiftLimit0 = 0; 
allSizex0 = zeros(1, 2) + sizex; 
allSizey0 = zeros(1, 2) + sizey; 
numImage0 = 1; 
SUM1map0 = cell(1, numOrient); 
MAX1map0 = cell(1, numOrient); 
tracking0 = cell(1, numOrient); 
deformedTemplate0{1} = single(zeros(sizex, sizey));
for (orient = 1 : numOrient)
       SUM1map0{1, orient} = single(zeros(sizex, sizey));
       MAX1map0{1, orient} = single(zeros(sizexSubsample, sizeySubsample));
       tracking0{1, orient} = single(zeros(sizexSubsample, sizeySubsample));
end
for (this = 1:numImage)
disp(['get example ' num2str(this)]);
for (orient = 1 : numOrient)
       Ccopy(SUM1map0{1, orient}, SUM1map{this, orient}, 0, 0, sizex, sizey, sizex, sizey); 
end
ClocalMaximumPooling(numImage0, allSizex0, allSizey0, numOrient,  ...
                     locationShiftLimit0, orientShiftLimit0, subsample, saturation, ...
                     SUM1map0, MAX1map0, tracking0);  
Clearn(numOrient, locationShiftLimit0, orientShiftLimit0, subsample, ... % about active basis  
       numElement, numImage0, sizex, sizey, ...
       SUM1map0, MAX1map0, tracking0, ... % about training images 
       halfFilterSize, Correlation, allSymbol(1, :), ... % about filters
       numStoredPoint, storedlambda, storedExpectation, storedLogZ, ... % about exponential model 
       selectedOrient, selectedx, selectedy, selectedlambda, selectedLogZ, ... % learned parameters
       commonTemplate, deformedTemplate0); % learned templates 
%% Iteration part 1: search nearest neighbors by template matching
for (it = 1 : numIteration)
close all; 
Cmatch(numOrient, subsample, ... % about active basis  
      numElement, numImage, sizex, sizey, MAX1map, ... % about training images 
      selectedOrient, selectedx, selectedy, selectedlambda, selectedLogZ, ...% learned parameters
      SUM2score) % output template matching score
%% Iteration part 2: find neighbors using current template
[SUM2scoreSort ind] = sort(SUM2score, 'descend'); 
for (k = 1 : numNeighbor)
     neighbork = ind(k); 
     for (orient = 1 : numOrient)
       Ccopy(SUM1mapLearn{k, orient}, SUM1map{neighbork, orient}, 0, 0, sizex, sizey, sizex, sizey); 
       Ccopy(MAX1mapLearn{k, orient}, MAX1map{neighbork, orient}, 0, 0, sizexSubsample, sizeySubsample, sizexSubsample, sizeySubsample);
       Ccopy(trackingLearn{k, orient}, tracking{neighbork, orient}, 0, 0, sizexSubsample, sizeySubsample, sizexSubsample, sizeySubsample); 
     end
end
Clearn(numOrient, locationShiftLimit, orientShiftLimit, subsample, ... % about active basis  
       numElement, numNeighbor, sizex, sizey, ...
       SUM1mapLearn, MAX1mapLearn, trackingLearn, ... % about training images 
       halfFilterSize, Correlation, allSymbol(1, :), ... % about filters
       numStoredPoint, storedlambda, storedExpectation, storedLogZ, ... % about exponential model 
       selectedOrient, selectedx, selectedy, selectedlambda, selectedLogZ, ... % learned parameters
       commonTemplate, deformedTemplate); % learned templates 
if (it==numIteration)
    meanSUM2score(this, 1) = mean(SUM2scoreSort(1:numNeighbor)); 
    allSUM2score(this, :) = SUM2score; 
    neighbor(this, :) = ind(1:numNeighbor); 
    allTemplate{this} = commonTemplate+0.; 
    Outthis; 
end
clear SUM2scoreSort; 
clear ind; 
end
save 'locallyLearnedTemplate' allTemplate allSUM2score meanSUM2score sizex sizey numImage;
end
OutAllTemplate; 


