
% Adapted by Anthony Ripa from code written by Prof Alex Berg

% debug level
debug_level = 3;

% window size
ws = 25;

sskip = 1;

% num_locs_per_image
num_locs_per_image = 100;

% num_locs_per_image_hard
num_locs_per_image_hard = 200;

% max_non_faces
max_non_faces = 10000;




%%%
%%%  load in example faces
%%%
%D = dir('example_faces/*.png');
face_dir = 'faces';
D = dir([face_dir '/*.jpg']);
D = D(1:55);

example_faces = zeros(ws,ws,length(D));

fc=0;  % initialize face count
for i=1:length(D),
    try 
        I = imread([[face_dir '/'],D(i).name]);
        I = im2double(I);  % represent elements as doubles instead of uint8
        I = rgb2gray(I);   % convert color to gray scale [rows x cols x 3] to [ rows x cols ]
        if (size(I,1)~=ws)||(size(I,2)~=ws),
            I = imresize(I,[ws,ws]);  % resize to a standard size
        end
        I = max(0,min(1,I)); %just fixing any values
                             %imresize might have produced outside [0,1]
        fc = fc + 1;
        example_faces(:,:,fc)=I;  % stack up the images in a 3d matrix
        if (debug_level > 2),  % optionally show the images as you read them in 
            figure(1);
            hold off
            imagesc(I);
            colormap(gray)
            % pause(0.2);
        end
    catch
    end
end

fprintf('\nloaded %d faces.\n',fc);

k = mean(example_faces,3)
imagesc(k);
title('mean face');
imwrite(k,'mean.jpg')


%%%
%%%  load in example nonfaces
%%%
D = dir('example_nonfaces/*.jpg');


example_non_faces = zeros(ws,ws,max_non_faces);

nfc = 0;  % number of non-faces
for i=1:length(D),
    I = imread(['example_nonfaces/',D(i).name]);
    I = im2double(I);  % represent elements as doubles instead of uint8
    I = rgb2gray(I);   % convert color to gray scale [rows x cols x 3] to [ rows x cols ]
    I = max(0,min(1,I)); %just fixing any values
                         %imresize might have produced outside [0,1]
    % pick random locations
    [h,w]=size(I);
    r = round(0.5+rand(num_locs_per_image)*(h-ws));
    r = max(1,min(h-ws+1,r));
    c = round(0.5+rand(num_locs_per_image)*(w-ws));
    c = max(1,min(w-ws+1,c));
    for j=1:num_locs_per_image,
        if (nfc==max_non_faces),
            break
        end
        nfc= nfc+1;
        example_non_faces(:,:,nfc) = I(r:r+ws-1,c:c+ws-1);
        if (debug_level > 5),  % optionally show the images as you read them in 
            figure(1);
            hold off
            imagesc(example_non_faces(:,:,nfc));
            colormap(gray)
            drawnow;  % force redraw
        end
    end
    if (nfc>max_non_faces),
        break
    end
end


fprintf('loaded %d faces and %d non-faces.\n',fc,nfc);

if (nfc>max_non_faces),
    fprintf('Exceeded max non-faces\n');
end




% interface with libsvm 
% http://www.csie.ntu.edu.tw/~cjlin/libsvm/
%
% install it: 
%   unpack, in root of installation type make.  
%   then in the matlab directory, from inside of matlab type make
%
%
% we need to add the path here:

% this path will depend on where you put libsvm
addpath ../libsvm-3.11/matlab

face_rp = randperm(fc);
perm_example_faces = example_faces(:,:,face_rp);

non_face_rp = randperm(nfc);
perm_example_non_faces = example_non_faces(:,:,non_face_rp);

% hold out some example for validation
n_face_train= fc-0;
n_non_face_train = nfc-100;

train = cat(3,perm_example_faces(:,:,1:n_face_train),perm_example_non_faces(:,:,1:n_non_face_train));
train = reshape(train,[ws*ws size(train,3)])';

train_labels = cat(1,ones(n_face_train,1),zeros(n_non_face_train,1));


 model = svmtrain(train_labels,train,'-s 0 -t 2 -c 1 -g .2'); %rbf
% model = svmtrain(train_labels,train,'-s 0 -t 0 -c 1'); %linear

% NonFace Accuracy = 1 - Alpha
fprintf('1 - Alpha = NonFace ');
test = cat(3,perm_example_non_faces(:,:,n_non_face_train+1:end)); 
test = reshape(test,[ws*ws size(test,3)])';
test_labels = cat(1,zeros(nfc-n_non_face_train,1));
[acc,aaa,bbb]=svmpredict(test_labels,test,model);

% Face Accuracy = 1 - Beta = Power
fprintf('Power = 1 - Beta = Face ');
test = cat(3,perm_example_faces(:,:,n_face_train+1:end)); 
test = reshape(test,[ws*ws size(test,3)])';
test_labels = cat(1,ones(fc-n_face_train,1));
[acc,aaa,bbb]=svmpredict(test_labels,test,model);

% Combined Accuracy
fprintf('Mixed ');
test = cat(3,perm_example_faces(:,:,n_face_train+1:end),perm_example_non_faces(:,:,n_non_face_train+1:end)); 
test = reshape(test,[ws*ws size(test,3)])';
test_labels = cat(1,ones(fc-n_face_train,1),zeros(nfc-n_non_face_train,1));
[acc,aaa,bbb]=svmpredict(test_labels,test,model);


% look at the support vectors
figure(5);
montage(permute(reshape(full(model.SVs'),[ws,ws,size(model.SVs,1)]),[1 2 4 3]));
colormap(gray)
title('ordered support vectors for initial model');
% the corresponding weights are in model.sv_coef


%%%
%%%  Try a test image again
%%%



I = imread('test2.jpg');
I = im2double(I);  % represent elements as doubles instead of uint8
I = rgb2gray(I);   % convert color to gray scale [rows x cols x 3]
                   % to [ rows x cols ]

if (max(size(I))>1000),
    I = imresize(I,0.5);
end

scales = sqrt(2).^[-3:3];
scales = [.8 .8];

for sc_index=1:length(scales),  % loop over scales
    Is = imresize(I,scales(sc_index));
    Is = max(0,min(1,Is));
    [h,w]=size(Is);
    % figure out offsets for sliding window locations, every sskip pixels
    h_samples = 1:sskip:h-ws+1;
    w_samples = 1:sskip:w-ws+1;
    nh_samples = length(h_samples);
    nw_samples = length(w_samples);
    n = floor(nh_samples*nw_samples);
    crops=zeros(ws*ws,n);
    pos = zeros(2,n);
    count =0;
    %extract the patches
    for i=h_samples,
        for j=w_samples,
            count = count + 1;
            c = Is(i:i+ws-1,j:j+ws-1);
            crops(:,count) = c(:);
            pos(:,count)=[i,j]';
        end
    end
    % evaluate the model on the extracted patches
    [cl,acc,clv]=svmpredict(zeros(count,1),crops(:,1:count)',model);
    % display the results for this scale 
    figure(4);
    clf
    subplot(2,1,1);
    result_im = reshape(clv,[nw_samples nh_samples])';
    imagesc(result_im);
    colormap(gray)
    subplot(2,1,2);
    imagesc(Is);
    imwrite(Is,'Is.png')
    hold on
    [sv,si]=sort(-clv);
    sv = -sv;
    for i=1:10,
        plot(pos(2,si(i)),pos(1,si(i)),'g*');
        plot(pos(2,si(i))+ws,pos(1,si(i))+ws,'g*');
        plot(pos(2,si(i))+ws,pos(1,si(i)),'g*');
        plot(pos(2,si(i)),pos(1,si(i))+ws,'g*');
    end
    drawnow
    sv(1:10)'
    all_Is{sc_index}  = Is;
    all_scale{sc_index} = scales(sc_index);
    all_result_im{sc_index} = result_im;
    all_clv{sc_index} = clv;
    all_pos{sc_index} = pos;
end

figure(3)
imagesc(all_result_im{2}); axis image; colormap(gray); colorbar
imwrite(all_result_im{2},'all.jpg')


