%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                 HOMEWORK #3 - EE148                                %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% Analyze images with both PCA and Fischer
clear
directories = {'Faces','Motorbikes'}; nd = size(directories,2);
SCALE = 1; %% Use SCALE=1/4 to subsample images by a factor of 4

%% Select subregion where the image actually contain the object
%% You can (and should) modify those values.
good_rows = [0.15 0.95];
good_cols = [0.2 0.8];

%% Proportion of the data allocated to training and testing
TRAIN = 4/5;
TEST = 1-TRAIN;

%% Read images from files and place them into matrix A
i=0; %% index into data matrix
A = []; Train = []; Test = [];
for d=1:nd,
    dir_name = directories{d};
    dir_cont = dir(dir_name); %% Contents of directory
    for f=3:length(dir_cont),
        file_name = dir_cont(f).name;
        if length(findstr(file_name,'jpg'))>0,
            i=i+1;
            ima = imread([dir_name '/' file_name],'jpg');
            [rr,cc] = size(ima);
            ima = ima(ceil(rr*good_rows(1)):ceil(rr*good_rows(2)),ceil(cc*good_cols(1)):ceil(cc*good_cols(2))); %% Select good region of the image
            ima = imresize(ima,SCALE,'bilinear'); %% Subsample the image
	    
		%% 1. Insert here code for grayscale normalization if so desired
            
	    figure(100); imagesc(ima); colormap(gray);
            A(:,i) = double(ima(:));
            Y(i) = d; %% Label of data (assume that data in a directory all belong to the same class)
        end;
    end;

    %% Keep track of how many examples are in each class and decide which images to use for training and which for testing
    if d>1, N(d) = i-N(d-1); else N(d) = i; end; %% Number of examples in each group
    I{d} = find(Y==d); %% index of each group
    i_max = round(TRAIN*N(d)); i_perm = randperm(N(d));
    Train_c{d} = I{d}(i_perm(1:i_max)); Train = [ Train Train_c{d}];
    Test_c{d} = I{d}(i_perm((i_max+1):N(d))); Test = [ Test Test_c{d} ];
end;


%% Calculate and display mean images on the whole dataset, just for the fun of seeing how they look
[rows,cols] = size(ima); %% Assume that all images have same size
mean_face = reshape(mean(A(:,I{2})')',rows,cols);
mean_motorbike = reshape(mean(A(:,I{1})')',rows,cols);
mean_face = reshape(mean(A')',rows,cols);
figure(1); 
subplot(1,3,1); imagesc(mean_face); colormap(gray); title('Mean face face');
subplot(1,3,2); imagesc(mean_motorbike); colormap(gray); title('Mean motorbike face');
subplot(1,3,3); imagesc(mean_face); colormap(gray); title('Mean face');


%%
%% 2. See if can separate faces and motorbikes with SVD
%%
mean_face_train = reshape(mean( A(:,Train)' )',rows,cols);
A0 = A(:,Train) - mean_face_train(:)*ones(1,length(Train)); %% Subtract mean

%% The following 2 lines of code compute the SVD of the matrix A0
%% Prove that they are equivalent to [U,S,V] = svd(A); 
%% Guess why we take this complicated route
[U1,S,V] = svd(A0' * A0);
U = A0 * V1 * inv(S1);

%% Display the first 16 principal components
figure(2); semilogy(diag(S),'o');
figure(333); 
for ef=1:16,
    subplot(4,4,ef); imagesc(reshape(U(:,ef),rows,cols)); colormap(gray); axis off; title(num2str(ef));
end;

dims = [8 12]; %% Pick two dimensions... try your favorite ones
mean_face_test = reshape(mean( A(:,Test)' )',rows,cols);
A0_test_faces = A(:,Test_c{2})-mean_face_test(:)*ones(1,length(Test_c{2}));
A0_test_motorbikes = A(:,Test_c{1})-mean_face_test(:)*ones(1,length(Test_c{1}));
A0_train_faces = A(:,Train_c{2})-mean_face_train(:)*ones(1,length(Train_c{2}));
A0_train_motorbikes = A(:,Train_c{1})-mean_face_train(:)*ones(1,length(Train_c{1}));
sub_faces = U(:,dims)' * A0_test_faces; 
sub_motorbikes = U(:,dims)' * A0_test_motorbikes;
figure(4); plot(sub_faces(1,:),sub_faces(2,:),'b+',sub_motorbikes(1,:),sub_motorbikes(2,:),'ro');
title('Motorbikes (red o) and faces (blue +)'); 
xlabel(['PC ' num2str(dims(1))]); ylabel(['PC ' num2str(dims(2))]);

%% 3. Here you need to add code to estimate class-conditional densities for faces and motorbikes
%% and plots the decision boundary for equal priors for faces and motorbikes.
%% Also code for plotting the ROC curve as you vary the threshold (or the priors).
%% You may pick the best
%% threshold either assuming that the prior probabilities are equal, or you may use the ROC to
%% select your favorite. You may use more than 2 dimensions if you wish.

%% See if can separate Faces from Motorbikes using mean face and motorbike faces
mean_face_train = mean(A(:,length(Train_c{1})+[1:length(Train_c{2})])')';
mean_motorbike_train = mean(A(:,1:length(Train_c{1}))')';
sub_faces = [mean_face_train mean_motorbike_train]' * A0_test_faces; 
sub_motorbikes = [mean_face_train mean_motorbike_train]' * A0_test_motorbikes;
figure(5); clf; plot(sub_faces(1,:),sub_faces(2,:),'b+',sub_motorbikes(1,:),sub_motorbikes(2,:),'ro');
title('Motorbikes (red o) and faces (blue +)'); 
xlabel('Mean faces'); ylabel('Mean motorbikes');

%% Again show class-conditional densities (with ellipses), the decision boundary, the ROC




%% 4. Now try Fisher. The code available on the homework page (fisherLD.m) will not work in this case... why?
%% Fix the code so that it works (the fix is very elementary, but understanding what to do is
%% not so easy... think about the meaning of the rank of a matrix and try and understand what
%% is the rank of the svd.
%% Once you have fixed the code Fisher will still not work well (excellent performance on training,
%% terrible performance on testing, as discussed in class). Why? What can you do to make
%% things better?
%II = [I{1} I{2}];
LDs  = fisherLD_PP ( A0, II(Train));
%LDs  = fisherLD ( A0,[ones(1,length(Train_c{1})), 2*ones(1,length(Train_c{2}))] );
figure(6);
for fld=1:16,
    subplot(4,4,fld); imagesc(reshape(LDs(:,fld),rows,cols)); colormap(gray); axis off; title(num2str(fld));
end;

dims = [1 2]; 
LD_faces = LDs(:,dims)'  * A0_test_faces;   
LD_motorbikes = LDs(:,dims)' * A0_test_motorbikes; 
figure(7); subplot(1,2,1); plot(LD_faces(1,:),LD_faces(2,:),'b+',LD_motorbikes(1,:),LD_motorbikes(2,:),'ro');
title('Test set: Motorbikes (red o) and faces (blue +)'); 
xlabel(['LD ' num2str(dims(1))]); ylabel(['LD ' num2str(dims(2))]);


LD_faces = LDs(:,dims)'  * A0_train_faces;   
LD_motorbikes = LDs(:,dims)' * A0_train_motorbikes;
figure(7); subplot(1,2,2); plot(LD_faces(1,:),LD_faces(2,:),'b+',LD_motorbikes(1,:),LD_motorbikes(2,:),'ro');
title('Training set: Motorbikes (red o) and faces (blue +)'); 
xlabel(['LD ' num2str(dims(1))]); ylabel(['LD ' num2str(dims(2))]);

%% Again show all the above. 


%%% 5. Try your own solution here!!!