% INF5300, spring 2016 - Lab - Linear feature transforms - Code related to the PCA part

load ~inf5300/www_docs/data/tm.mat

% Place the images in one Nx6 matrix
X = double([Converted_from_raw_data_1(:) Converted_from_raw_data_2(:) Converted_from_raw_data_3(:) Converted_from_raw_data_4(:) Converted_from_raw_data_5(:) Converted_from_raw_data_6(:)]);

N = size(X,1); % Let N be the number of samples/pixels

% Subtract the mean
X = X-repmat(mean(X), size(X,1), 1);

% Create sample covariance matrix / scatter matrix
R = (X'*X)/N;

% Get the eigenvectors and eigenvalues of R (note the ascending ordering of the eigenvalues we get from eig())
[A, D] = eig(R);

% Use the columns of eigenvectors as linear-transform weights to get the PCA components
Y = X*A;

% Create sample covariance matrix / scatter matrix of our transformed variables
Ry = (Y'*Y)/N

% Reshape the data into a six-banded 2D image
imPca = reshape(Y, size(Converted_from_raw_data_1,1), size(Converted_from_raw_data_1,2), 6);

% Have a look at the "pca" images
figure;
for i=1:6
    subplot(3,2,i);
    imagesc(imPca(:,:,i));
    title(sprintf('lambda = %.2f', D(i,i)));
end
colormap gray

% Study the weights behind the principle component image
figure;
stem(A(:,end))
title('Feature-weights for the first principle component');
xlabel('Original feature number')
ylabel('Weight')


% See if our simple loop gives us our principle eigenvector:
r = randn(6,1);
for i=1:100
    r = R*r;
    r = r/norm(r);
end
r        % Remember that we do not differentiate between r and -r
A(:,end) % The eigenvector corresponding to the highest eigenvalue

% Remove the r component from R and do procedure again:
R2 = R - D(end,end)*r*r';
r = randn(6,1);
for i=1:100
    r = R2*r;
    r = r/norm(r);
end
r
A(:,end-1) % The eigenvector corresponding to the next-to-highest eigenvalue


