% 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 
numElement = 45; % number of Gabors in active basis   
numIteration = 5;  % number of iterations
numMixture = 3; 
mkdir('working'); 
%% 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, localOrNot, localHalfx, localHalfy, thresholdFactor); % 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(numImage, numOrient); 
for (img = 1:numImage)
    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, numImage); % templates for training images 
for (img = 1 : numImage)
    deformedTemplate{img} = single(zeros(sizex, sizey));  
end
%% 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
for (cate = 1:numMixture)
   this = floor(rand(1)*numImage) + 1; 
   disp(['start learning from example ' num2str(this)]); tic
   for (orient = 1 : numOrient)
       Ccopy(SUM1map0{1, orient}, SUM1map{this, orient}, sizex, sizey); 
   end
   deformedTemplate0{1} = single(zeros(sizex, sizey));
   CsharedSketch(numOrient, locationShiftLimit0, orientShiftLimit0, subsample, ... % about active basis  
       numElement, numImage0, sizex, sizey, SUM1map0, ... % 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 
   figure; imshow(-commonTemplate, []); 
   templateName = ['working/learnedTemplate' num2str(cate)]; 
   save(templateName, 'selectedOrient', 'selectedx', 'selectedy', 'selectedlambda', 'selectedLogZ'); 
end
%% Iteration part 1: search nearest neighbors by template matching
for (it = 1 : numIteration)
SUM2score = zeros(numImage, numMixture);
for (cate = 1:numMixture)
    templateName = ['working/learnedTemplate' num2str(cate)]; 
    load(templateName); 
    SUM2scorecate = zeros(numImage, 1); 
    CgetSUM2(numOrient, subsample, ... % about active basis  
      numElement, numImage, sizex, sizey, MAX1map, ... % about training images 
      selectedOrient, selectedx, selectedy, selectedlambda, selectedLogZ, ...% learned parameters
      SUM2scorecate) % output template matching score
    SUM2score(:, cate) = SUM2scorecate+0.; 
end
num = zeros(numMixture, 1); 
example = zeros(numMixture, numImage); 
for (img = 1 : numImage)
     [max_SUM2, cate_img] = max(SUM2score(img, :));
     num(cate_img) = num(cate_img) + 1; 
     example(cate_img, num(cate_img)) = img; 
end
num
for (cate = 1:numMixture)
  for (k = 1:num(cate))
    img = example(cate, k); 
    for (orient = 1 : numOrient)
       Ccopy(SUM1mapLearn{k, orient}, SUM1map{img, orient}, sizex, sizey); 
    end
  end
%% Iteration part 2: re-learn the template from nearest neighbors
CsharedSketch(numOrient, locationShiftLimit, orientShiftLimit, subsample, ... % about active basis  
       numElement, num(cate), sizex, sizey, SUM1mapLearn(1:num(cate), :), ... % 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 
   templateName = ['working/learnedTemplate' num2str(cate)]; 
   save(templateName, 'selectedOrient', 'selectedx', 'selectedy', 'selectedlambda', 'selectedLogZ'); 
if (it==numIteration)
bx =0; by = 0; cc = 3;
col = 6; row = floor(num(cate)/col)+1;
Iout = zeros(row*sizex+(row-1)*(2*bx+cc), 2*col*sizey+(col-1)*(2*by+cc))+255; 
towrite = -commonTemplate; 
towrite = (255 * (towrite-min(towrite(:)))/(max(towrite(:))-min(towrite(:))));
Iout(1:sizex, 1:sizey) = towrite;
k = -1;  
for (r = 1 : row)
    for (c = 1 : col)
        k = k+1; 
        if ((k>=1)&&(k<=num(cate)))
            img = example(cate, k); 
            towrite = Ifill{img}; 
            towrite = (255 * (towrite-min(towrite(:)))/(max(towrite(:))-min(towrite(:))));
            Iout((r-1)*sizex+1+(r-1)*(2*bx+cc): r*sizex+(r-1)*(2*bx+cc), (c-1)*2*sizey+1+(c-1)*(2*by+cc): (c-1)*2*sizey+sizey+(c-1)*(2*by+cc)) = towrite; 
            towrite =  -deformedTemplate{k}; 
            towrite = (255 * (towrite-min(towrite(:)))/(max(towrite(:))-min(towrite(:))));
            Iout((r-1)*sizex+1+(r-1)*(2*bx+cc) : r*sizex+(r-1)*(2*bx+cc), (c-1)*2*sizey+sizey+(c-1)*(2*by+cc)+1 : (c-1)*2*sizey+sizey+(c-1)*(2*by+cc)+sizey) =towrite;  
        end
    end
end
for (r = 1 : row-1)
  Iout(r*sizex+(r-1)*(2*bx+cc)+bx+1:r*sizex+(r-1)*(2*bx+cc)+bx+cc, :) = 0; 
end
for (c = 1 : col-1)
  Iout(:, c*(sizey*2)+(c-1)*(2*by+cc)+by+1:c*(sizey*2)+(c-1)*(2*by+cc)+by+cc) = 0; 
end
  outepsgif(Iout, ['mixture_' num2str(cate)]); 
end %i if 
end % cate
end % it





