%--  3:32 PM 9/11/08 --%
%
% Matlab -- PCA ANALYSIS
%

% download from
% web.media.mit.edu/~maov/classes/vision08/assignments/ps3/freiburg_train.mat
load freiburg_train

% dtensor_train is a 5-way array
% 75 people were photographed from 6 view points under 
% 6 illumination conditions.
% Each image has 80 rows and 107 cols
size(dtensor_train)
% let's look at a few images
figure; imshow(reshape(dtensor_train(1,1,1,:,:),80,107),[])
figure; imshow(reshape(dtensor_train(1,3,1,:,:),80,107),[])

% Each image will be vectorized using the function reshape
D=reshape(dtensor_train,75,6,6,80*107);
% The array will be turned into a matrix that has 80*107 pixels and 2700 images
X=reshape(D,75*6*6,80*107)';
size(X)

% Let's display a couple of images that are currently in matrix X
x=reshape(X(:,1),80,107);
figure;imshow(x,[])
x=reshape(X(:,2),80,107);
figure;imshow(x,[])

% mean image
mnd=mean(X,2);
D=double(X)-repmat(mnd,1,2700);
imshow(reshape(uint8(mnd),80,107),[])

% two ways of computing the covariance
cv=cov(D);
cv=D*D'/(2700-1);

% Compute the eigenvectors.
% These can be compute using the function svd on the covariance of D,
% or on matrix D. It is comp. more efficient to do an svd(D,0) than an svd(cv,0).
% Check the size of cv and size of D.

[U,S,V]=svd(D,0);
% Display some of the eigenvectors
figure; imshow(reshape(U(:,1),80,107),[])
figure; imshow(reshape(U(:,2),80,107),[])

% Notice that the eigenvectors displayed above are orthonormal to each other
% (Two ways of computing the dot product) 
dot(U(:,1),U(:,2))
(U(:,1)'*U(:,2))

% Let's display the top 5 eigenvectors
allaxis=[]; 
for v=1:5,
  mn=min(U(:,v));
  mx=max(U(:,v));
  im=(U(:,v)-mn)/(mx-mn); %img normalization values btw 0&1
  im=reshape(im,80,107);
  allaxis=[allaxis im];
end
figure; imshow(allaxis,[])

% Computing the coefficient matrix
C=U'*D;

% New Representation for image 1
c1=C(:,1);

% Image 1 reconstruction
d1=U*c1;
figure;imshow(uint8(reshape(d1,80,107)),[])
figure;imshow(uint8(reshape(d1+mnd,80,107)),[])

% Reconstruction after dimensionality reduction
% Coefficient c1 reduced from 2700 -> 1000
d1=U(:,1:1000)*c1(1:1000,1);
figure; imshow(uint8(reshape(d1+mnd,80,107)),[])

% Coefficient c1 reduced from 2700 -> 100
d1=U(:,1:100)*c1(1:100,1);
imshow(uint8(reshape(d1+mnd,80,107)),[])
diff=D(:,1)-d1;
error=sqrt(diff'*diff)
