% Lezione di lab.TTR  -21'30/03/2016- PCA e eigenfaces - Marco Cristani

clear all;
close all;


% % -------------------------------------------------------------
% % Eigenfaces
% % -------------------------------------------------------------

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 image; colormap gray;
%         drawnow;
%         pause(0.01);
    TMP(:,i)= tmp(:);
%         mont(:,:,1,i)= tmp;
end

% figure;
% montage(mont);
% set(gcf,'name','Gallery');

%A lezione abbiamo visto una slide che riporta delle raccomandazioni su
%come debbano essere organizzate le immagini per poter applicare in
%sicurezza eigenfaces. Dando un'occhiata alle raccomandazioni, quali di
%esse sono soddisfatte dall'ORL? E quali no?


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

[U,lambda]              =   eigen_training(A);

% Andare a leggere e capire il codice di eigen_training, facendo mapping
% con la lezione vista in classe.

figure;
subplot(211)
plot(lambda); title('Autovalori')
subplot(212)
plot(cumsum(lambda)/sum(lambda)); title('Informazione modellata')


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

% 
% % SOL 2) Effettuare prove di proiezione con T autovettori

% -1D  T=1-
% Prendere tutte le immagini a disposizione, e proiettarle su una
% dimensione, avendo poi cura di dar loro e visualizzare un'etichetta 
% testuale che ne identifichi la persona (vi ricordo, abbiamo 40 persone,
% ognuna di esse ha 10 istanze)

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

figure; set(gcf,'name','Visualizzazione 1D dei volti');
scatter(omega(1,:),ones(1,M)); hold on

l           =   reshape(repmat([1:40],10,1),400,1);
for i=1:M
    text(omega(1,i),1+0.01,num2str(l(i)))
end

% Verificare che il primo autovettore è quello lungo cui i punti
% si disperdono maggiormente. A questo proposito, provare a fare 
% proiezioni 1D su autovettori diversi dal primo.

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

figure; set(gcf,'name','Visualizzazione 1D dei volti su autovettore minore');
scatter(omega(1,:),ones(1,M)); hold on

l           =   reshape(repmat([1:40],10,1),400,1);
for i=1:M
    text(omega(1,i),1+0.01,num2str(l(i)))
end



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

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

l           =   reshape(repmat([1:40],10,1),400,1);
% for i=1:M
%     text(omega(1,i)+25,omega(2,i)+25,num2str(l(i)))
% end

for i=1:M
    text(omega(1,i)+25,omega(2,i)+25,num2str(i))
end



%% 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;

figure; set(gcf,'name','Similarita nello spazio 2D = spazio orig?');
subplot(1,3,1); imagesc(reshape(TMP(:,313),r,c));axis image
title('313 simile a 123')
subplot(1,3,2); imagesc(reshape(TMP(:,123),r,c));axis image; 
title('123')
subplot(1,3,3); imagesc(reshape(TMP(:,97),r,c));axis image; colormap gray;
title('dissimile da 313 e 123')


%% Analisi massima somiglianza e diversita' nel dataset;
% Valutare le proiezioni massimamente vicine e lontane, osservarne le
% immagini corrispondenti, variando la dimensionalità della proiezione

T           =   150;
omega       =   U(:,1:T)'*A; % proiezione;
for i=1:M
    for j=1:M
        dista(i,j)= norm(omega(:,i)'-omega(:,j)');
    end
end

figure; imagesc(dista)
colorbar
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 50D = 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
hhh=figure; 
rec = [];
rico = 123;
omega = [];
for t=1:400
    T              =   t;
    omegatmp       =   U(:,T)'*A(:,rico); % proiezione;
% %%Annullamento di una componente: cosa accade?
%     if t==4
%         omegatmp = -omegatmp;
%     end
    omega          =   [omega; omegatmp];
    
    rec = media + U(:,1:T)*omega;
    subplot(1,3,1); imagesc(reshape(TMP(:,rico),r,c)); colormap gray
    title('originale'); axis image
     subplot(1,3,2); imagesc(reshape(rec,r,c)); colormap gray
    title(strcat('ricostruita con ',num2str(t),'-autovettori')); axis image
    subplot(1,3,3); imagesc(reshape(U(:,t),r,c)); colormap gray
    title(strcat('eigenface n:',num2str(t),'|a_i=',num2str(omega(t)))); axis image
    drawnow; pause (0.001); 
    mov(t) = getframe(hhh)
end







%% Analisi delle singole basi

figure;
for i=1:25
    subplot(5,5,i); imagesc(reshape(U(:,i),[r,c]));
    colormap gray; axis image; axis off;
%     title(num2str(lambda(i),'%0.5g'));
    title(lambda(i));
end


Tan           =   [1]; %base analizzata
omega       =   U(:,Tan)'*A; % proiezione su un'unica base;


[thrash, ind] = sort(omega,'ascend');

figure; set(gcf,'name','Analisi singole basi');
conta = 0;

for i=1:1:25
    conta =conta +1;
subplot(5,5,conta); 
imagesc(reshape(TMP(:,ind(i)),r,c));axis image; 
title(num2str(ind(i)))
end
colormap gray




%% Prove di analisi di robustezza della ricostruzione



% Aggiungo ad un'immagine del rumore impulsivo
clear U
clear lambda
clear A
ind_img               =   110
good                    =   setdiff([1:M], ind_img);
TMPnew                  =   TMP;
TMPnew(:,ind_img)     =   [];
media                   =   mean(TMPnew,2);
A(:,:)                  =   TMPnew-repmat(media,1,M-1);

[U,lambda]              =   eigen_training(A);

img_noise = TMP(:,ind_img);
figure; set(gcf,'name','Prove di robustezza al rumore')
subplot(131)
imagesc(reshape(img_noise,r,c)); axis image; colormap gray; axis off
title('img originale')
img_noise = imnoise(uint8(img_noise),'salt & pepper', 0.02);
subplot(132)
imagesc(reshape(img_noise,r,c)); axis image; colormap gray; axis off
title('img con noise')

omega_noise = U'*(double(img_noise)-media);
rec = media + U*omega_noise;
subplot(133)
imagesc(reshape(rec,r,c)); axis image; colormap gray; axis off
title('ricostruzione');

omega       =   U'*A; % per avere le proiezioni della dimensionalita'
                      %  di partenza
distanza_orig = norm(omega_noise - omega(:,100))

%Prova di occlusione;
img_occ = TMP(:,80);
figure; set(gcf,'name','Prove di robustezza ad occlusione')
subplot(131)
imagesc(reshape(img_occ,r,c)); axis image; colormap gray; 
title('img originale')

img_occ = reshape(img_occ,r,c);
img_occ(50:70,20:40)=10;
img_occ(50:70,55:75)=10;
subplot(132)
imagesc(reshape(img_occ,r,c)); axis image; colormap gray; 
title('img occlusa');

omega_occ = U'*(double(img_occ(:))-media);
rec = media + U*omega_occ;
subplot(133)
imagesc(reshape(rec,r,c)); axis image; colormap gray; axis off
title('ricostruzione');


% Prendere un'immagine e proiettarla su 2D con PCA, 
% ottenendo w. Aggiungere del rumore all'immagine, 
% proiettarla, ottenendo w2. Come sono relazionati tra di 
% loro w e w2 in termini di distanza?

ind_img   =   110
img_noise = TMP(:,ind_img);
figure; set(gcf,'name','Prove di robustezza al rumore')
subplot(131)
imagesc(reshape(img_noise,r,c)); axis image; colormap gray; axis off
title('img originale')
img_noise = imnoise(uint8(img_noise),'salt & pepper', 0.02);
subplot(132)
imagesc(reshape(img_noise,r,c)); axis image; colormap gray; axis off
title('img con noise')

img_occ = TMP(:,ind_img)
img_occ = reshape(img_occ,r,c);
img_occ(50:70,20:40)=10;
img_occ(50:70,55:75)=10;

subplot(133)
imagesc(reshape(img_occ,r,c)); axis image; colormap gray; 
title('img occlusa');

omega = U'*(double(TMP(:,ind_img))-media);
omega_noise = U'*(double(img_noise)-media);
omegaall =  U'*A;

figure; set(gcf,'name','Visualizzazione 2D dei volti');
scatter(omegaall(1,:),omegaall(2,:));hold on
scatter(omega(1,:),omega(2,:),'r');
scatter(omega_noise(1,:),omega_noise(2,:),'k');
omega_occ = U'*(double(img_occ(:))-media);
scatter(omega_occ(1,:),omega_occ(2,:),'g');
%% Riconoscimento tramite eigenfaces


clear all
close 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));
    TMP(:,i)= tmp(:);
end

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






%Eseguo C cicli
C=50;
hrec = figure; set(gcf,'name','RICONOSCIMENTO VOLTI');
label = reshape(repmat([1:40],10,1),400,1);
conf = zeros(40,40);
for iter=1:400
    % Estrazione probe (elemento incognito)
    index_probe = iter;%round(unifrnd(1,M)); %indice probe
    probe = TMP(:,index_probe);
    subplot(211); imagesc(reshape(probe,r,c)); title('probe');...
        colormap gray, axis image;
    newTMP = TMP;
    newTMP(:,index_probe) = [];
    newlabel = label;
    newlabel(index_probe) = [];
    % 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           =   399;
    omega_g       =   U(:,1:T)'*newA; % proiezione della gallery;
    omega_p     =   U(:,1:T)'*(probe-media);
    num_neigh = 5;
    class = knnclassify(omega_p',omega_g',newlabel,num_neigh);
    figure(hrec); subplot(212); imagesc(reshape(TMP(:,class*10-1),r,c));
    axis image, colormap gray;
    accuracy(iter)=class==label(index_probe);
     title(strcat('matched?:',num2str(accuracy(iter))))
    fprintf('run %i: accuracy = %i\n',iter,accuracy(iter));
    conf(label(index_probe),class)=conf(label(index_probe),class)+1;
    drawnow
    pause(0.01);
end

%mostrare accuratezza e matrice di confusione del classificatore
%soprastante


sum(accuracy)/50

% %% RICOSTRUZIONE
% % Estrazione probe (elemento incognito)
% num_tolti =10;
% 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           =   300;
% 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(recon(:,i),r,c)); title('rec ');colormap gray, axis image;
% end