% Joe Jupin
% Homework 7
%modified by Longin Jan Latecki, latecki@temple.edu
%October 2004

%  clear all variables
clear all;
%close all;

noPCcomp=10; %number of PC components


% read images
% names of images
fpics = {'carpet.bmp',
        'fur.bmp',
        'grass.bmp',
        'gravel.bmp',
        'knit.bmp',
        'tartan.bmp'};
fpics = fpics';
% get size of gpics
fnum = numel(fpics);
% read images into double matrix
for i=1:fnum
    f(:,:,i) = imread(char(fpics(i)));
end
% get gpics dimensions
[ipixs jpixs kpics] = size(f);

% obtain a stack of 24 subsize x subsize images from original images
subsize = 32;
for i=1:fnum
    f1(:,:,(i*4-1)) = f(1:subsize,1:subsize, i);
    f1(:,:,(i*4-2)) = f((ipixs-subsize+1):ipixs,1:subsize, i);
    f1(:,:,(i*4-3)) = f(1:subsize,(jpixs-subsize+1):jpixs, i);
    f1(:,:,(i*4)) = f((ipixs-subsize+1):ipixs,(jpixs-subsize+1):jpixs, i);
end
% obtain number and dimensions of images to be processed
[f1ipixs f1jpixs f1kpics] = size(f1);
Q = f1kpics;

for i=1:f1kpics
a=f1(:,:,i);
% X(i,:)=a(:);
[x, x_pos, imagerows, imagecols] = cut_evzr5(a, subsize/2, subsize/2);
X(i,:)=x';
end
X=double(X);
meanx=mean(X);
meanX=repmat(meanx,size(X,1),1);
normalX=X-meanX;

[pc, latent, explained] = pcacov(cov(normalX'));%transpose trick
%pc=normalX'*pc;%transpose trick

PCAproj=pc(:,1:noPCcomp);


projX=(PCAproj'*normalX'); %projection of the images

recX=(PCAproj*projX)'+ meanX; %reconstruction of the images

% for i=1:Q %reformating reconstructed images
%     pv = recX(i,:);
%     pv = reshape(pv, subsize, subsize);
%     recXim(:,:,i) = pv;
% end
% 
% 
% 
% projEV=(PCAproj'*PCAproj); %projection of the eigenvectors
% 
% recEV=(PCAproj*projEV)'; %reconstruction of the eigenvectors
% 
% for i=1:noPCcomp %reformating reconstructed eigenvectors images
%     pev = recEV(i,:);
%     pev = reshape(pev, subsize, subsize);
%     recEVim(:,:,i) = pev;
% end



% show original images
figure;
for i=1:f1kpics
    subplot(4,Q/4,i);
    imshow(f1(:,:,i));
end

% % show reconstructed images
% figure;
% for i=1:Q
%     subplot(4,Q/4,i);
%     imshow(recXim(:,:,i),[]);
% end
% 
% % show eigenimages
% figure;
% for i=1:noPCcomp
%     subplot(4,round(noPCcomp/4),i);
%     imshow(recEVim(:,:,i),[]);
% end


%section on quallity of the reconstruction

%distance computation
ds=recX;
    s=size(ds);
    ds_distances = zeros(s(1),s(1));
        for i=1:s(2)
            a=ds(:,i)*ds(:,i)';
            b=repmat(ds(:,i),1,s(1)).^2;
            c=repmat(ds(:,i)',s(1),1).^2;
            ds_distances = b+c-2*a + ds_distances;
        end
        ds_dist=sqrt(ds_distances);

figure; imagesc(ds_dist)


%we compute the retrieval rate for each of the 6 classes 
%sumres contains a sumetric matrix with pariwise distances
symres=ds_dist;
for i=0:5
    for j=1:4
        [a,b]=sort(symres(i*4+j,:));
        c=b(1,1:4);
        q=find(i*4<c & c<i*4+5);
        retrat4(i*4+j)=size(q,2)/4;
    end
    retrat=mean(retrat4);
end

retrat