% LEARNING ACTIVE BASIS AND ADABOOST TEMPLATE
%% 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 ClocalNormalizeDouble.c; % local normalization of type double
mex ClocalNormalize.c; % local normalization
mex Clearn.c; % learning by shared sketch algorithm
mex CgetMAX1.c; % local maximum pooling for MAX1 maps
mex Csqrt.c; % local maximum pooling for MAX1 maps
mex CdrawTemplate.c; % local normalization
epsilon = .1; % allowed correlation between selected Gabors 
subsample = 1; % 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 = 'ABvsADwomanMoreNeg';
%% create output folders 
if (exist([outputFolder 'eps'])==0)
    mkdir([outputFolder 'eps']); 
else
    delete([outputFolder 'eps/*.*']); 
end
if (exist([outputFolder 'png'])==0)
    mkdir([outputFolder 'png']); 
else 
    delete([outputFolder 'png/*.*']); 
end
if (exist([outputFolder '.html']))
    delete([outputFolder '.html']); 
end
fid = fopen([outputFolder '.html'], 'wt'); 
%% Load in training images
imageOption = 1; resizeFactor = .3; % 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; 
%% Compute SUM1 maps by Gabor filtering
disp(['start filtering']); tic
if (imageOption == 3) 
    [SUM1map, Ifill] = ApplyFilterfftFill(I, allFilter, localOrNot, localHalfx, localHalfy, thresholdFactor); % SUM1 maps by Gabor filtering
else 
    SUM1map = ApplyFilterfftSame(I, allFilter, localOrNot, localHalfx, localHalfy, -1, thresholdFactor);% SUM1 maps by Gabor filtering
end
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
selectedPolar = zeros(1, numElement);
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']);
%% Display active basis templates 
fprintf(fid, '%s\n', ['<a href="' outputFolder 'Code.zip">Code and data</a> <br>']);
outMessage =  ['==> Number of training images = ' num2str(numImage) ...
                          '; Number of elements = ' num2str(numElement) ...
                          '; Length of Gabor = ' num2str(halfFilterSize*2+1) ' pixels' ...
                          '; Local normalization or not = ' num2str(localOrNot) ...
                          '; Range of displacement = ' num2str(locationShiftLimit) ' pixels' ...
                          '; Subsample rate = ' num2str(subsample) ' pixel' ...
                          '; Image height and width = '  num2str(sizex) ' by ' num2str(sizey) ' pixels ' ...
                          '<br>'];
fprintf(fid, '%s\n', outMessage);
disp(outMessage);    
fprintf(fid, '%s\n', ['<hr>']);
heightstr = '"height=170> ';  
for (numSketch = 5:5:numElement)
    CdrawTemplate(numOrient, numSketch, sizex, sizey, subsample, halfFilterSize, allSymbol(1, :), selectedOrient, selectedx, selectedy, selectedPolar, commonTemplate); 
    pic = -double(commonTemplate); showImage; 
    saveas(gcf, [outputFolder 'eps/AB' num2str(1000+numSketch) '.eps'], 'eps');
    saveas(gcf, [outputFolder 'png/AB' num2str(1000+numSketch) '.png'], 'png');
    fprintf(fid, '%s\n', ['<IMG SRC="' outputFolder 'png/AB' num2str(1000+numSketch)  '.png' heightstr]);
end
fprintf(fid, '%s\n', ['<br>']);
fprintf(fid, '%s', ['Number of sketches in the above active basis templates are, respectively']);
for (numSketch = 5:5:numElement)
    fprintf(fid, '%s', [', ' num2str(numSketch)]);
end
fprintf(fid, '%s\n', ['<br>']);
%% Filtering and maxing positive images 
disp(['start filtering positive images and computing MAX1 maps']); tic
if (imageOption == 3) 
    [positiveSUM1map, positiveIfill] = ApplyFilterfftFill(I, allFilter, localOrNot, localHalfx, localHalfy, thresholdFactor); % SUM1 maps by Gabor filtering
else 
    positiveSUM1map = ApplyFilterfftSame(I, allFilter, localOrNot, localHalfx, localHalfy, -1, thresholdFactor);% SUM1 maps by Gabor filtering
end
positiveNumImage = numImage; 
allSizex = zeros(1, positiveNumImage)+sizex; 
allSizey = zeros(1, positiveNumImage)+sizey; 
Csqrt(positiveNumImage, allSizex, allSizey, numOrient, positiveSUM1map);
sizexSubsample = floor(sizex/subsample); 
sizeySubsample = floor(sizey/subsample); 
positiveMAX1map = cell(positiveNumImage, numOrient);
for (img = 1 : positiveNumImage)
    for (orient = 1:numOrient) 
        positiveMAX1map{img, orient} = single(zeros(sizexSubsample, sizeySubsample));
    end
end
CgetMAX1(positiveNumImage, allSizex, allSizey, numOrient, locationShiftLimit, orientShiftLimit, subsample, ...
         positiveSUM1map, positiveMAX1map); 
disp(['filtering time: ' num2str(toc) ' seconds']);
%% Crop negative images for adaboost training
maxSize = 800; boundary = 20; % resize the original image to maxSize and remove the boundary 
naturalImageFolder = 'naturalImage'; % folder of training images  
samplingRate = .3;  % sampling the image patches with probability samplingRate
downScaleFactor = 1/2.; % repeatedly resize the image for cropping at multiple resolution 
lowerVar = .01; % skip flat image patches of small variance 
nonOverlapFactor = 1/2.; % gap between consecutive image patches to be sampled and cropped 
increasex = floor(sizex*nonOverlapFactor); increasey = floor(sizey*nonOverlapFactor); 
CropNegativeImage; % crop negative image patches of (sizex, sizey) from natural images 
%% Filtering and maxing negative images 
disp(['start filtering negative images and computing MAX1 maps']); tic
negativeSUM1map = ApplyFilterfftSame(negativeI, allFilter, localOrNot, localHalfx, localHalfy, -1, thresholdFactor);% SUM1 maps by Gabor filtering
allSizex = zeros(1, negativeNumImage)+sizex; 
allSizey = zeros(1, negativeNumImage)+sizey; 
Csqrt(negativeNumImage, allSizex, allSizey, numOrient, negativeSUM1map);
negativeMAX1map = cell(negativeNumImage, numOrient);
for (img = 1:negativeNumImage)
    for (orient = 1:numOrient) 
        negativeMAX1map{img, orient} = single(zeros(sizexSubsample, sizeySubsample));
    end
end
CgetMAX1(negativeNumImage, allSizex, allSizey, numOrient, locationShiftLimit, orientShiftLimit, subsample, ...
         negativeSUM1map, negativeMAX1map); 
disp(['filtering time: ' num2str(toc) ' seconds']);
%% Train adaboost template 
dataWeight = single(0.5 * [ones(1, positiveNumImage)/positiveNumImage, ones(1, negativeNumImage)/negativeNumImage]); % starting from balanced weights for adaboost 
adnumElement = 2*numElement; % maximum number of selected weak classifiers 
numGridPoint = numStoredPoint; % number of grid points to search for threshold, to be the same as number of grid in lambda function of active basis 

positiveFeatureOnly = 1; 
if (positiveFeatureOnly>0.)
   mex adaboostPosOnly.c;         
   [selectedXs selectedYs selectedOs selectedThresholds selectedPolars lambdas]...
        = adaboostPosOnly(positiveMAX1map, negativeMAX1map,...
                           dataWeight,numGridPoint,halfFilterSize,1,adnumElement);
   ADname = 'ADposonly'; 
   announce = 'Features include only positive ones. ';
else
   mex adaboostLearning.c;         
   [selectedXs selectedYs selectedOs selectedThresholds selectedPolars lambdas]...
        = adaboostLearning(positiveMAX1map, negativeMAX1map,...
                           dataWeight,numGridPoint,halfFilterSize,1,adnumElement);
   ADname = 'AD'; 
   announce = 'Features include both positive and negative ones. ';
end
selectedPolarsAD = (double(selectedPolars)-.5)*2; 

%% Display active basis templates 
for (numSketch = 5:5:adnumElement)
    CdrawTemplate(numOrient, numSketch, sizex, sizey, subsample, halfFilterSize, allSymbol(1, :), double(selectedOs), double(selectedXs), double(selectedYs), double(selectedPolarsAD), commonTemplate); 
    pic = -double(commonTemplate); showImage; 
    saveas(gcf, [outputFolder 'eps/' ADname num2str(1000+numSketch) '.eps'], 'eps');
    saveas(gcf, [outputFolder 'png/' ADname num2str(1000+numSketch) '.png'], 'png');
    fprintf(fid, '%s\n', ['<IMG SRC="' outputFolder 'png/' ADname num2str(1000+numSketch)  '.png' heightstr]);
end
fprintf(fid, '%s\n', ['<br>']);
fprintf(fid, '%s', [announce 'Number of negative images is ' num2str(negativeNumImage) '. Dark sketches correspond to positive features. Light sketches correspond to negative features. Number of sketches in the above adaboost templates are, respectively']);
for (numSketch = 5:5:adnumElement)
    fprintf(fid, '%s', [', ' num2str(numSketch)]);
end
fprintf(fid, '%s\n', ['<br>']);
