%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                            %
% Practical assignment: Face detection                                       %
%                                                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                            %
%   Part 1:  load and prepare image the data                                 %
%   Part 2:  learn linear and RBF SVM face classifier                        %
%   Part 3:  scanning window face detection with non-max suppression         %
%                                                                            %
%   Note: the code below relies on libSVM library (included). If you         %         
%   encounter problems with running 'svmtrain' and 'svmpredict' functions,   %
%   go to directory ../../practical-face-detection/libsvm-mat-3.0-1 and      %
%   run 'make' either inside Matlab (Windows) or shell (Linux).              %
%                                                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


basepath='../../practical-face-detection';
imgpath=[basepath '/images'];
respath=[basepath '/results'];

addsvmpath
colormap(gray)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% Part 1: Loading and normalizing training images %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%% load training images
%%%%%%%%%%%%%%%%%
posfname=[imgpath '/possamples.mat'];
negfname=[imgpath '/negsamples.mat'];
load(posfname,'possamples'); npos=size(possamples,3);
load(negfname,'negsamples'); nneg=size(negsamples,3);  
possamples=double(possamples);
negsamples=double(negsamples);
fprintf('load  %6d positive samples from %s\n',npos,posfname)
fprintf('load  %6d negative samples from %s\n\n',nneg,negfname)


%%%%%%%%%%%%%%%% positive and negative training images are now 
%%%%%%%%%%%%%%%% inside 'possamples' and 'possamples' 3D arrays
%%%%%%%%%%%%%%%% of the size 24x24xN where
%%%%%%%%%%%%%%%%  - 24x24 is the pixel size of croped faces
%%%%%%%%%%%%%%%%  - N is the number of samples
%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%% It's always good to look at the data
%%%%%%%%%%%%%%%% let's display a few positive and negative samples
%%%%%%%%%%%%%%%%
possampleimage=[];
negsampleimage=[];
for i=1:10
  indpos=randperm(npos);
  indneg=randperm(nneg);
  posfaces=reshape(shiftdim(possamples(:,:,indpos(1:10)),1),24*10,24);
  negfaces=reshape(shiftdim(negsamples(:,:,indneg(1:10)),1),24*10,24);
  possampleimage=cat(2,possampleimage,posfaces);
  negsampleimage=cat(2,negsampleimage,negfaces);
end

clf
fprintf('Showing images of random Positive samples\n')
showimage(possampleimage'); title('positive samples')
fprintf('press a key...'), pause, fprintf('\n')
fprintf('Showing images of random Negative samples\n')
showimage(negsampleimage'); title('negative samples')
fprintf('press a key...'), pause, fprintf('\n')



%%%%%%%%%%%%%%%% we can also display an average face and non-face images
%%%%%%%%%%%%%%%%
clf, 
fprintf('Showing average images for Positive and Negative samples\n')
subplot(1,2,1), title('average positive image')
showimage(mean(possamples(:,:,1:2000),3)); 
subplot(1,2,2), title('average negative image')
showimage(mean(negsamples(:,:,1:2000),3)); 
fprintf('press a key...'), pause, fprintf('\n')


%%%%%%%%%%%%%%%% Positive training images have been manually aligned
%%%%%%%%%%%%%%%% by x-y coordinates. However, they may have large
%%%%%%%%%%%%%%%% variation in the amplitude. To remove this variation
%%%%%%%%%%%%%%%% we normalize pixel values of all samples to zero mean
%%%%%%%%%%%%%%%% and variance value=1
%%%%%%%%%%%%%%%%
possamples=meanvarpatchnorm(possamples);
negsamples=meanvarpatchnorm(negsamples);


%%%%%%%%%%%%%%% SVM training and classification assumes each sample
%%%%%%%%%%%%%%% in the form of a vector. Let's flatten training images
%%%%%%%%%%%%%%% into one vector per sample
%%%%%%%%%%%%%%%
xsz=size(possamples,2);
ysz=size(possamples,1);
Xpos=transpose(reshape(possamples,ysz*xsz,npos));
Xneg=transpose(reshape(negsamples,ysz*xsz,nneg));


%%%%%%%%%%%%%%% make vectore with sample labels:
%%%%%%%%%%%%%%% +1 for positives
%%%%%%%%%%%%%%% -1 for negatives
%%%%%%%%%%%%%%%
ypos=ones(npos,1);
yneg=-ones(nneg,1);
  
  
%%%%%%%%%%%%%%% separate data into the training and validation set
%%%%%%%%%%%%%%% 
ntrainpos=1000;
ntrainneg=1000;
indpostrain=1:ntrainpos; indposval=indpostrain+ntrainpos;
indnegtrain=1:ntrainneg; indnegval=indnegtrain+ntrainneg;

Xtrain=[Xpos(indpostrain,:); Xneg(indnegtrain,:)];
ytrain=[ypos(indpostrain); yneg(indnegtrain)];
Xval=[Xpos(indposval,:); Xneg(indnegval,:)];
yval=[ypos(indposval); yneg(indnegval)];

% free memory
clear possamples negsamples



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% Part 2: Train test SVM classifier, analyze parameters%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%% set SVM parameters: linear kernel with C-parameter specified
%%%%%%%%%%%%%%%

%-t kernel_type : set type of kernel function (default 2)
%	0 -- linear: u'*v
%	1 -- polynomial: (gamma*u'*v + coef0)^degree
%	2 -- radial basis function: exp(-gamma*|u-v|^2)
%	3 -- sigmoid: tanh(gamma*u'*v + coef0)
%-d degree : set degree in kernel function (default 3)
%-g gamma : set gamma in kernel function (default 1/num_features)
%-c cost : set the parameter C of C-SVC, epsilon-SVR, and nu-SVR (default 1)
%-v : verbose mode

C=0.001;
libsvmparamstr=sprintf('-t 0 -c %1.5f',C);


%%%%%%%%%%%%%%% train linear SVM
%%%%%%%%%%%%%%%
fprintf('Training SVM classifier with %d pos. and %d neg. samples...',sum(ytrain==1),sum(ytrain~=1))
model = svmtrain(ytrain,Xtrain,libsvmparamstr);
yalpha=model.sv_coef;
Xsup=model.SVs;
b=model.rho;
fprintf(' -> %d support vectors\n',size(Xsup,1))

 

%%%%%%%%%%%%%%% get prediction for training and validation samples
%%%%%%%%%%%%%%%
fprintf('Running evaluation... ')
[predtrain,acctrain,conftrain] = svmpredict(ytrain,Xtrain,model);
[predval,  accval,  confval]   = svmpredict(yval,  Xval,  model);
fprintf('Training accuracy: %1.3f; validation accuracy: %1.3f\n',acctrain(1),accval(1))
fprintf('press a key...'), pause, fprintf('\n')



%%%%%%%%%%%%%%% **************************************************************
%%%%%%%%%%%%%%% **************************************************************
%%%%%%%%%%%%%%% *                                                            *
%%%%%%%%%%%%%%% *                       EXERCISE 1:                          *
%%%%%%%%%%%%%%% *                                                            *
%%%%%%%%%%%%%%% *       Compute SVM hyper-plane W from support vectors       *
%%%%%%%%%%%%%%% *                                                            *
%%%%%%%%%%%%%%% **************************************************************
%%%%%%%%%%%%%%% **************************************************************
%%%%%%%%%%%%%%%                                                             
%%%%%%%%%%%%%%% Linear SVM classifier can be expressed in the form of a 
%%%%%%%%%%%%%%% hyper-plane W in the original feature space, i.e. 24*25=576-D
%%%%%%%%%%%%%%% space of face-image pixel values. Let's construct W.
%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%% Recall that W can be expressed in terms of support vectors
%%%%%%%%%%%%%%% and their coefficients as W = sum_i [alpha_i*y_i*X_i]
%%%%%%%%%%%%%%% The solution above provides support vecotors in 'Xsup'
%%%%%%%%%%%%%%% and alpha_i*y_i in 'yalpha'. 
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% TODO:
%%%%%%%%%%%%%%% 1.1 Compute W from Xsup and yalpha
%%%%%%%%%%%%%%% 1.2 re-compute confidence for training and validation samples values
%%%%%%%%%%%%%%%     as distances between sample vectors X and W using bias 'b' as 
%%%%%%%%%%%%%%%     conf=X*W-b
%%%%%%%%%%%%%%% 1.3 re-compute classification accuracy using true sample labels 'y'
%%%%%%%%%%%%%%%     as acc=mean((conf>0)*2-1==y) for training and validation samples


% 1.1
W=transpose(sum(Xsup.*repmat(yalpha,1,size(Xsup,2))));

% 1.2 
conftrainnew = Xtrain*W-b;
confvalnew   = Xval*W-b;

% 1.3
acctrainnew = mean((conftrainnew>0)*2-1==ytrain);
accvalnew   = mean((confvalnew>0)*2-1==yval);
fprintf('Training and validation accuracy re-estimated from W,b: %1.3f; %1.3f\n',acctrainnew,accvalnew)
fprintf('press a key...'), pause, fprintf('\n')

%%%%%%%%%%%%%%% **************************************************************
%%%%%%%%%%%%%%% ******************* END of EXERCISE 1 ************************
%%%%%%%%%%%%%%% **************************************************************



%%%%%%%%%%%%%%% The values of the hyper-plane W are
%%%%%%%%%%%%%%% the weights of individual pixel values and can be displayed
%%%%%%%%%%%%%%% as an image. Let's construct W from support vectors
%%%%%%%%%%%%%%%
clf, showimage(reshape(W,24,24))





%%%%%%%%%%%%%%% **************************************************************
%%%%%%%%%%%%%%% **************************************************************
%%%%%%%%%%%%%%% *                                                            *
%%%%%%%%%%%%%%% *                       EXERCISE 2:                          *
%%%%%%%%%%%%%%% *                                                            *
%%%%%%%%%%%%%%% *        Optimize C parameter on the validation set          *
%%%%%%%%%%%%%%% *                                                            *
%%%%%%%%%%%%%%% **************************************************************
%%%%%%%%%%%%%%% **************************************************************
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% SVM generalization ability depends on the constant C
%%%%%%%%%%%%%%% A standard practice for choosing C is to empirically
%%%%%%%%%%%%%%% select the best value of C by maximizing performance on
%%%%%%%%%%%%%%% the validation set. (Note that C or any other parameters
%%%%%%%%%%%%%%% should never be optimized on the final test set.)
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% TODO:
%%%%%%%%%%%%%%% 2.1 train linear SVM as above for different values of C
%%%%%%%%%%%%%%%     and select the best model 'modelbest' maximizing
%%%%%%%%%%%%%%%     accuracy on the validation set. Hint: try exponentially
%%%%%%%%%%%%%%%     distribution values, e.g. [1000 100 10 1 .1 .01 .001 .0001 .00001];
%%%%%%%%%%%%%%% 2.2 When training for different C values, compute W,b
%%%%%%%%%%%%%%%     and visualize W as an image (see above)
%%%%%%%%%%%%%%%
Call=[1000 100 10 1 .1 .01 .001 .0001 .00001];
accbest=-inf; 
modelbest=[];
for i=1:length(Call)
  C=Call(i);
  libsvmparamstr=sprintf('-t 0 -c %1.5f',C);
  fprintf('Training SVM classifier with C=%1.5f...',C)
  model = svmtrain(ytrain,Xtrain,libsvmparamstr);
  yalpha=model.sv_coef;
  Xsup=model.SVs;
  b=model.rho;
  W=transpose(sum(Xsup.*repmat(yalpha,1,size(Xsup,2))));
  fprintf(' -> %d support vectors\n',size(Xsup,1))
  
  W=transpose(sum(Xsup.*repmat(yalpha,1,size(Xsup,2))));
  conftrain = Xtrain*W-b;
  confval   = Xval*W-b;
  acctrain = mean((conftrain>0)*2-1==ytrain);
  accval   = mean((confval>0)*2-1==yval);
  fprintf('Training accuracy: %1.3f; validation accuracy: %1.3f\n',acctrain,accval)
  
  % select best model
  if accbest<accval
    accbest=accval;
    modelbest=model;
    Cbest=C;
  end
  
  clf, showimage(reshape(W,24,24))
  title(sprintf('C=%1.5f | Training accuracy: %1.3f; validation accuracy: %1.3f',C,acctrain,accval))
  drawnow
end
fprintf(' -> Best accuracy %1.3f for C=%1.5f\n',accbest,Cbest)
fprintf('press a key...'), pause, fprintf('\n')

%%%%%%%%%%%%%%% **************************************************************
%%%%%%%%%%%%%%% ******************* END of EXERCISE 2 ************************
%%%%%%%%%%%%%%% **************************************************************




%%%%%%%%%%%%%%% get W,b for the best model
%%%%%%%%%%%%%%%
yalphabest=modelbest.sv_coef;
Xsupbest=modelbest.SVs;
bbest=modelbest.rho;
Wbest=transpose(sum(Xsupbest.*repmat(yalphabest,1,size(Xsupbest,2))));

  

%%%%%%%%%%%%%%% W of the best classifier "looks like" and average face image
%%%%%%%%%%%%%%% Is it sufficient to construct W as the mean of of positive
%%%%%%%%%%%%%%% samples??? Let's try it
%%%%%%%%%%%%%%%


Walt=transpose(mean(Xtrain(find(ytrain==1),:),1));
clf, showimage(reshape(Walt,24,24)), drawnow

balt=0;
conftrainnew = Xtrain*Walt-balt;
confvalnew   = Xval*Walt-balt;
acctrainnew = mean((conftrainnew>0)*2-1==ytrain);
accvalnew   = mean((confvalnew>0)*2-1==yval);
fprintf('Training and validation accuracy estimated from W = "average image": %1.3f; %1.3f\n',acctrainnew,accvalnew)
fprintf('press a key...'), pause, fprintf('\n')


%%%%%%%%%%%%%%% train SVM with RBF kernel
%%%%%%%%%%%%%%%
C=10; 
gamma=1/500;
libsvmparamstr=sprintf('-t 2 -c %1.5f -g %1.5f',C,gamma);
fprintf('Training SVM with RBF kernel for %d pos. and %d neg. samples...',sum(ytrain==1),sum(ytrain~=1))
modelrbf = svmtrain(ytrain,Xtrain,libsvmparamstr);
fprintf(' -> %d support vectors\n',size(model.SVs,1));


%%%%%%%%%%%%%%% get prediction for training and validation samples
%%%%%%%%%%%%%%%
fprintf('Running evaluation... ')
[predtrain,acctrain,conftrain] = svmpredict(ytrain,Xtrain,modelrbf);
[predval,  accval,  confval]   = svmpredict(yval,  Xval,  modelrbf);
fprintf('RBF training accuracy: %1.3f; RBF validation accuracy: %1.3f\n',acctrain(1),accval(1))
fprintf('press a key...'), pause, fprintf('\n')



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% Part 3: Scanning window face detection %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%% load test image
%%%%%%%%%%%%%%%
imgfname='img3.jpg';
img=imread([imgpath '/' imgfname]);
showimage(img), title(imgfname)
gimg=mean(img,3);

%%%%%%%%%%%%%%% decompose the image into sub-windows
%%%%%%%%%%%%%%%
[ysz,xsz,csz]=size(img);
wxsz=24; wysz=24;
[x,y]=meshgrid(1:xsz-wxsz,1:ysz-wysz);
bbox=[x(:) y(:) x(:)+wxsz-1 y(:)+wysz-1];
n=size(bbox,1);
imgcropall=zeros(wysz,wxsz,length(x(:)));
for i=1:size(bbox,1)
  imgcropall(:,:,i)=cropbbox(gimg,bbox(i,:));
end



%%%%%%%%%%%%%%% Linear SVM classifier can be evaluated efficiently by
%%%%%%%%%%%%%%% using dot-product between image patches and W. Before
%%%%%%%%%%%%%%% computing the confidence we need to normalize and reshape
%%%%%%%%%%%%%%% test image patches.

imgcropall=meanvarpatchnorm(imgcropall);
X=transpose([reshape(imgcropall,wysz*wxsz,n)]);
conf=X*Wbest-bbest;


%%%%%%%%%%%%%%% display most confident detections
%%%%%%%%%%%%%%%
n=20;
[vs,is]=sort(conf,'descend');
clf, showimage(img), showbbox(bbox(is(1:n),:))
title(sprintf('%d best detections',n),'FontSize',14)
fprintf('press a key...'), pause, fprintf('\n')





%%%%%%%%%%%%%%% **************************************************************
%%%%%%%%%%%%%%% **************************************************************
%%%%%%%%%%%%%%% *                                                            *
%%%%%%%%%%%%%%% *                       EXERCISE 3:                          *
%%%%%%%%%%%%%%% *                                                            *
%%%%%%%%%%%%%%% *        Non-maxima suppression of multiple responses        *
%%%%%%%%%%%%%%% *                                                            *
%%%%%%%%%%%%%%% **************************************************************
%%%%%%%%%%%%%%% **************************************************************
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% Scanning-window style classification of image patches typically
%%%%%%%%%%%%%%% results in many multiple responses around the target object.
%%%%%%%%%%%%%%% A standard practice to deal with this is to remove any detector
%%%%%%%%%%%%%%% responses in the neighborhood of detections with locally maximal
%%%%%%%%%%%%%%% confidence scores (non-maxima suppression or NMS). NMS is
%%%%%%%%%%%%%%% usually applied to all detections in the image with confidence
%%%%%%%%%%%%%%% above certain threshold.
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% TODO:
%%%%%%%%%%%%%%% 3.1 Try out different threshold values to pre-selected windows
%%%%%%%%%%%%%%%     passed to the NMS stage, see parameter 'confthresh' below.
%%%%%%%%%%%%%%% 3.2 Try out different threshold values for NMS detections,
%%%%%%%%%%%%%%%     see parameter 'confthreshnms'
%%%%%%%%%%%%%%% 3.3 Try detection and with different thresholds for different
%%%%%%%%%%%%%%%     included images: 'img1.jpg', 'img2.jpg', 'img3.jpg', 'img4.jpg'
%%%%%%%%%%%%%%% 
confthresh=0;
indsel=find(conf>confthresh);
[nmsbbox,nmsconf]=prunebboxes(bbox(indsel,:),conf(indsel),0.2);


%%%%%%%%%%%%%%% display detections above threshold after non-max suppression
%%%%%%%%%%%%%%% 
confthreshnms=1;
clf, showimage(img)
indsel=find(nmsconf>confthreshnms);
showbbox(nmsbbox(indsel,:),[1 1 0],regexp(num2str(nmsconf(indsel)'),'\d+\.\d+','match'));
title(sprintf('%d NMS detections above threshold %1.3f',size(nmsbbox,1),confthreshnms),'FontSize',14)
fprintf('press a key...'), pause, fprintf('\n')


%%%%%%%%%%%%%%% **************************************************************
%%%%%%%%%%%%%%% ******************* END of EXERCISE 3 ************************
%%%%%%%%%%%%%%% **************************************************************



%%%%%%%%%%%%%%% Improve detection with non-linear RBF SVM classifier
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% Evaluation using RBF kernel takes much longer because we need to 
%%%%%%%%%%%%%%% compute distances to support vectors explicitly --> evaluate
%%%%%%%%%%%%%%% detection which have passed linear SVM model with low threshold
%%%%%%%%%%%%%%% 

confthresh=0;
indsel=find(conf>confthresh);
[predrbf,accrbf,confrbf] = svmpredict(ones(size(indsel)),X(indsel,:),modelrbf);
[nmsbboxrbf,nmsconfrbf]=prunebboxes(bbox(indsel,:),confrbf,0.2);


%%%%%%%%%%%%%%% display detections above threshold after non-max suppression
%%%%%%%%%%%%%%% 

confthreshnms=1;
clf, subplot(1,2,1), showimage(img)
indsel=find(nmsconf>confthreshnms);
showbbox(nmsbbox(indsel,:),[1 1 0],regexp(num2str(nmsconf(indsel)'),'\d+\.\d+','match'));
title('Linear','FontSize',14)

confthreshrbf=1;
subplot(1,2,2), showimage(img)
indsel=find(nmsconfrbf>confthreshrbf);
showbbox(nmsbboxrbf(indsel,:),[1 1 0],regexp(num2str(nmsconfrbf(indsel)'),'\d+\.\d+','match'));
title('RBF','FontSize',14)

