% Lezione di lab. 2 TTR  -11/04/2012- PCA e eigenfaces - Marco Cristani
%%% Simple PCA
clear all;
close all;

% -------------------------------------------------------------
% Esercizio 1.1; Riduzione di dimensionalita' di una gaussiana
% -------------------------------------------------------------


M1                                  =   [2,7]';
Cov1                                =   [0.2,0;0,2];
Data1                               =   mvnrnd(M1,Cov1,200);
h = figure; set(gcf,'name','PCA')
subplot(2,2,1); scatter(Data1(:,1),Data1(:,2)); axis([0, 12,0,12]);
title('Dati ibridi'); grid on;
[p(:,1),p(:,2)]                     =   cart2pol(Data1(:,1),Data1(:,2));
p(:,1)                              =   p(:,1)-pi/6;
[Data1(:,1),Data1(:,2)]             =   pol2cart(p(:,1),p(:,2));
subplot(2,2,2); scatter(Data1(:,1),Data1(:,2)); axis([0 12,0,12]);
title('Dati originali'); grid on;

Mu                                  =   mean(Data1)';
TMP                                 =   Data1'; % 
[num_dim, num_s]                    =   size(TMP);

A                                   =   TMP-repmat(Mu,1,num_s);
subplot(2,2,3); scatter(A(1,:)',A(2,:)'); title('Dati traslati'); 
 axis([-6 6,-6,6]);grid on;

S                                   =   A*A';
[E,lambda]                          =   eigs(S);

x =-10:10;
subplot(2,2,4); scatter(A(1,:)',A(2,:)'); hold on;
plot([E(1,1)*-10,E(1,1)*10],[E(2,1)*-10,E(2,1)*10],'r')
plot([E(1,2)*-10,E(1,2)*10],[E(2,2)*-10,E(2,2)*10],'b')
title('Eigenvectors'); axis equal; grid on;

%Controllate lambda....

% -------------------------------------------------------------
% Esercizio 1.2; Eigenfaces
% -------------------------------------------------------------
clear all;
dire    =    '.\images'; % ORL face database
list    =   dir(strcat(dire,'\*.bmp'));
M       =   size(list,1);
tmp     =   imread(strcat(dire,'\',list(1).name));
[r,c]   =   size(tmp);
TMP     =   zeros(r*c,M);
figure;
for i=1:M
    tmp = imread(strcat(dire,'\',list(i).name));
%     imagesc(tmp); axis equal; colormap gray;
%     drawnow;
%     pause(0.01);
    TMP(:,i)= tmp(:);
end

TMP                     =   single(TMP);
media                   =   mean(TMP,2);
A(:,:)                  =   TMP-repmat(media,1,M);

[U,lambda]              =   eigen_training(A);


%% Visualizzazione delle eigenfaces
figure;
for i=1:25
    subplot(5,5,i); imagesc(reshape(U(:,i),[r,c]));
    colormap gray; axis image; title(num2str(i));
end


% SOL 2) Effettuare prove di proiezione con T autovettori
% 2D
T           =   2;
omega       =   U(:,1:T)'*A; % proiezione;

figure; set(gcf,'name','Visualizzazione 2D dei volti');
scatter(omega(1,:),omega(2,:));

%% Analisi di similarita' tra immagini nello spazio ridotto:
%   Sono sufficienti 2 dimensioni per inferire somiglianze nello 
%   spazio originale? Guardare gli indici nello spazio ridotto, selezionare
%   dele proiezioni vicine e lontane, osservare visualmente i dati
%   originali corrispondenti;
hold on;
for i=1:M
    a=double(omega(1,i));
    b=double(omega(2,i));
    text(a,b,num2str(i));
end
figure; set(gcf,'name','Similarita nello spazio 2D = spazio orig?');
subplot(1,3,1); imagesc(reshape(TMP(:,331),r,c));axis image
title('331 simile a 135')
subplot(1,3,2); imagesc(reshape(TMP(:,135),r,c));axis image; 
title('135')
subplot(1,3,3); imagesc(reshape(TMP(:,127),r,c));axis image; colormap gray;
title('dissimile da 135 e 331')


%% Analisi massima somiglianza e diversita' nel dataset;
% Valutare le proiezioni massimamente vicine e lontane, osservarne le
% immagini corrispondenti.

for i=1:M
    for j=1:M
        dista(i,j)= norm(omega(:,i)'-omega(:,j)');
    end
end
dista = dista + diag(ones(1,M)*inf);
[distmin,Imin]= min(dista(:));
[sim_i,sim_j]=ind2sub([M,M],Imin);
[distmax,Imax]= max(dista(dista~=inf));
[dsim_i,dsim_j]=ind2sub([M,M],Imax);

figure; set(gcf,'name','Similarita nello spazio 2D = spazio orig?');
subplot(2,2,1); imagesc(reshape(TMP(:,sim_i),r,c));axis image
title(strcat(num2str(sim_i),'-simile-'))
subplot(2,2,2); imagesc(reshape(TMP(:,sim_j),r,c));axis image
title(strcat(num2str(sim_j),'-simile-'))
subplot(2,2,3); imagesc(reshape(TMP(:,dsim_i),r,c));axis image
title(strcat(num2str(dsim_i),'-dissimile-'))
subplot(2,2,4); imagesc(reshape(TMP(:,dsim_j),r,c));axis image
title(strcat(num2str(dsim_j),'-dissimile-'))
colormap gray

%% Analisi similarita' 1D; 
% Come prima, ma 1D

T           =   1;
omega       =   U(:,1:T)'*A; % proiezione;

for i=1:M
    for j=1:M
        dista(i,j)= norm(omega(:,i)'-omega(:,j)');
    end
end
dista = dista + diag(ones(1,M)*inf);
[distmin,Imin]= min(dista(:));
[sim_i,sim_j]=ind2sub([M,M],Imin);
[distmax,Imax]= max(dista(dista~=inf));
[dsim_i,dsim_j]=ind2sub([M,M],Imax);

figure; set(gcf,'name','Similarita nello spazio 1D = spazio orig?');
subplot(2,2,1); imagesc(reshape(TMP(:,sim_i),r,c));axis image
title(strcat(num2str(sim_i),'-simile-'))
subplot(2,2,2); imagesc(reshape(TMP(:,sim_j),r,c));axis image
title(strcat(num2str(sim_j),'-simile-'))
subplot(2,2,3); imagesc(reshape(TMP(:,dsim_i),r,c));axis image
title(strcat(num2str(dsim_i),'-dissimile-'))
subplot(2,2,4); imagesc(reshape(TMP(:,dsim_j),r,c));axis image
title(strcat(num2str(dsim_j),'-dissimile-'))
colormap gray

%% Analisi similarita' multi
% COme prima, 15 dimensioni
T           =   15;
omega       =   U(:,1:T)'*A; % proiezione;

for i=1:M
    for j=1:M
        dista(i,j)= norm(omega(:,i)'-omega(:,j)');
    end
end
dista = dista + diag(ones(1,M)*inf);
[distmin,Imin]= min(dista(:));
[sim_i,sim_j]=ind2sub([M,M],Imin);
[distmax,Imax]= max(dista(dista~=inf));
[dsim_i,dsim_j]=ind2sub([M,M],Imax);

figure; set(gcf,'name','Similarita nello spazio multi = spazio orig?');
subplot(2,2,1); imagesc(reshape(TMP(:,sim_i),r,c));axis image
title(strcat(num2str(sim_i),'-simile-'))
subplot(2,2,2); imagesc(reshape(TMP(:,sim_j),r,c));axis image
title(strcat(num2str(sim_j),'-simile-'))
subplot(2,2,3); imagesc(reshape(TMP(:,dsim_i),r,c));axis image
title(strcat(num2str(dsim_i),'-dissimile-'))
subplot(2,2,4); imagesc(reshape(TMP(:,dsim_j),r,c));axis image
title(strcat(num2str(dsim_j),'-dissimile-'))
colormap gray


%% Effettuare prove di ricostruzione, al variare del numero di eigenfaces
% figure; 
% for t=1:400
%     T           =   t;
%     omega       =   U(:,1:T)'*A; % proiezione;
%     rec = repmat(media,1,M) + U(:,1:T)*omega;
%     subplot(2,1,2); imagesc(reshape(rec(:,100),r,c)); colormap gray
%     title(strcat('ricostruita con ',num2str(t),'-autovettori')); axis image
%     subplot(2,1,1); imagesc(reshape(TMP(:,100),r,c)); colormap gray
%     title('originale'); axis image
%     pause(0.25)
% end

%% Riconoscimento tramite eigenfaces
% Estrazione probe (elemento incognito)
index_probe = round(unifrnd(1,M)); %indice probe
probe = TMP(:,index_probe);
hrec = figure; set(gcf,'name','RICONOSCIMENTO VOLTI');
subplot(211); imagesc(reshape(probe,r,c)); title('probe');colormap gray, axis image;
label = reshape(repmat([1:40],10,1),400,1);
newTMP = [TMP(:,1:index_probe-1),TMP(:,index_probe+1:end)]; %problema se probe e' l'ultimo elto;
newlabel = [label(1:index_probe-1);label(index_probe+1:end)];
% Creazione eigenvalues
media                   =   mean(newTMP,2);
newA                    =   newTMP-repmat(media,1,M-1);
[U,lambda]              =   eigen_training(newA);
% Proiezione incognito, decido il numero di compo.
T           =   150;
omega_g       =   U(:,1:T)'*newA; % proiezione della gallery;
omega_p     =   U(:,1:T)'*(probe-media);
class = knnclassify(omega_p',omega_g',newlabel,20);
figure(hrec); subplot(212); imagesc(reshape(TMP(:,class*10),r,c));
axis image, colormap gray; title('vote:')

%% RICOSTRUZIONE
% % Estrazione probe (elemento incognito)
% num_tolti =4;
% index_class= round(unifrnd(1,40)); %indice probe sulla classe
% probe = TMP(:,index_class*10-num_tolti+1:index_class*10);
% hrec = figure; set(gcf,'name','RICOSTRUZIONE VOLTI');
% for i=1:num_tolti
% subplot(2,num_tolti,i); imagesc(reshape(probe(:,i),r,c)); title('probe');colormap gray, axis image;
% end
% newTMP = [TMP(:,1:index_class*10-num_tolti),TMP(:,index_class*10+1:end)]; %prob se probe e' l'ultimo elto;
% % Creazione eigenvalues
% media                   =   mean(newTMP,2);
% newA                    =   newTMP-repmat(media,1,M-num_tolti);
% [U,lambda]              =   eigen_training(newA);
% % Proiezione incognito, decido il numero di compo.
% T           =   200;
% omega_p     =   U(:,1:T)'*(probe-repmat(media,1,num_tolti));
% recon       =   U(:,1:T)*omega_p +repmat(media,1,num_tolti);
% for i=1:num_tolti
% subplot(2,num_tolti,i+num_tolti); imagesc(reshape(probe(:,i),r,c)); title('rec ');colormap gray, axis image;
% end