% Always close everything beforehand
close all;
clear all;

%
% 1. Load all images
% Load in a sample image set (here we use lighting set as the example) with m images, and flatten them so the pixels in each image is a column vector, for instance using ImV=Im(:);. The image as a vector will be a quite long vector, e.g. q = 10,000 for a 100 x 100 image. Then put all the m column vectors into an q x m matrix M.
% The default "loadlightimages" script, along with the suggested indices for light images "1" - (29 37 24 50 11 61 2 64 5 56 17 45 33) are used.
%

% You can vary the number and X for different datasets
number = 1; 
X = [29 37 24 50 11 61 2 64 5 56 17 45 33];
num_frames = size(X, 2);
imsize = [480, 640];

% ########################################
% WRITE YOUR CODE HERE TO CONSTRUCT matrix A
% ########################################
M = zeros(prod(imsize),num_frames);

% ########################################

for i = 1:size(X, 2)
    idx = X(i);
    im = loadlightimages(number, idx, [imsize(2), imsize(1)])';

    % ########################################
    % WRITE YOUR CODE HERE TO flatten im 
    
    % ########################################
end

%
% 2. Perform a Normalized Eigen-vector/Principle Components Analysis (PCA) on your image set
% (Linear Subspace for a Set of Images) Formulate a linear system so that BX = Y, where B is the basis vector, Y contains the K eigenvectors, X is linearly spaced values for time representation.
% 2.0 Choose the No. Basis (K) to Save
% Expect this number to vary.

K = 10; % K

% 2.1 Mean Image
% Find the mean of all input image vectors and make Mz a matrix where each column is an image minus the mean image (meanIm). This can be achieved by the built in Matlab function mean, but pay attention to which dimension it takes the mean over. We will synthesize/interpolate the intensity variation of a pixel over time, so the mean should be over the different image instances of the same pixel, not the mean over the pixel intensities in one image.
% Hint: Use an outer product of a mean column vector and a row vector of ones to form a matrix of the same size as Mz with the mean image repeated in each column. Optionally, visualize the acquired mean image and denote its association to the image dataset. Now you can directly subtract this mean matrix from M.

% ########################################
% WRITE YOUR CODE HERE TO calculate a mean image Imean 
Imean = mean(M,2);

% ########################################
imshow(reshape(Imean,imsize));

% ########################################
% WRITE YOUR CODE HERE TO form Mz
Mz = M-Imean*ones(1,num_frames);

% ########################################

%
% 2.2 Normalized Covariance Matrix C
% Calculate the covariance matrix C = (Mz' * Mz) / m;
% where m is the number of images in our dataset.
% Note that C will be symmetric and positive definite, with size m x m.

% ########################################
% WRITE YOUR CODE HERE TO
% Calculate normalized covariance matrix C here


% ########################################

%
% 2.3 Perform PCA
% Get the first K eigenvectors Y of C using eigs in matlab. (K=6 will be needed below, but for better quality try K=10 for the light, and K=20 for the objects. )
% The eigenvectors we calculate could be used as a basis to reconstruct the columns of Mz'. The columns of Mz' each represent how one pixel in the image changes over time.
% Note that by default eigs will only give the eigenvectors corresponding to the 6 highest eigenvalues.
%

% ########################################
% WRITE YOUR CODE HERE TO
% Use eigs to acquire eigenvectors and eigenvalues. Store the eigenvectors
% in a matrix Y.


% ########################################

%
% 2.4 Calculate and save the basis images
% We will calculate a basis B that can be used to reconstruct entire images from our dataset.
% The eigenvectors Y calculated above can be used as a basis for the transpose of our image matrix Mz'. Therefore, we have
% {\mathbf{M}}_{\mathbf{z}}^{\mathbf{T}} ={\textrm{YB}}^{\mathbf{T}}
% for some matrix B'.
% We want a basis for Mz instead of Mz'. We can take the transpose of both sides of the equation above.
% {\mathbf{M}}_{\mathbf{z}} ={\textrm{BY}}^{\mathbf{T}}
% Y is an orthonormal matrix of eigenvectors, so {\left({\mathbf{Y}}^{\mathbf{T}} \right)}^{-1} =\mathbf{Y}. So we can calculate B as
% {\mathbf{M}}_{\mathbf{z}} \mathbf{Y}={\textrm{BY}}^{\mathbf{T}} \mathbf{Y}
% {\mathbf{M}}_{\mathbf{z}} \mathbf{Y}=\mathbf{B}
% Be sure to remember testing your own generated basis using script from the previous exercise.

% Calculate an image basis B 

B = Mz*Y;

% As above we have Mz = B*Y'. We can reconstruct our images from the basis B
% with the coefficicents in Y'. Let's save Y' as Y (this is the same Y as in ex4).

Y = Y';

% To stay consistent with the default saved X
X = [-120 -95 -70 -50 -25 -10 0 10 25 50 70 95 120];    % For light set

Imean = reshape(Imean, imsize);
save('light_pca_ex5.mat', "Y", "X", "B", "imsize", "Imean");

%
% 3. Visualization
% 3.1 Visualize your Basis Vectors
% Transform the resulting basis vectors (columns of B) into images of the same size as the input images, and display the first 5 (use absolute value). Make sure you get the basis vectors and images corresponding to the largest eigenvalues.
%

% ########################################
% WRITE YOUR CODE HERE TO
% Visualize your Image Basis


% ########################################

% 3.2 Plot the coefficients of Y
% Plot a graph of the corresponding coefficients (rows of Y); these are the functions we want to interpolate using polynomials (one poly for each row of Y).

% ########################################
% WRITE YOUR CODE HERE TO
% Visualize your coefficients


% ########################################

%
% 5. Normalized SVD
% A normalized version of SVD.
%

% ########################################
% WRITE YOUR CODE HERE TO
% Perform SVD over the normalized Anorm


% ########################################

%
% Calculate and save the basis images for SVD.

% Select the best num_basis
% ########################################
% WRITE YOUR CODE HERE TO
% How do we select K for basis B and Y?


% ########################################
save('light_svd_ex5.mat', "Y", "X", "B", "imsize", "Imean");

%
% Comparison: SVD vs. PCA
% Which method gives a better quality in image reconstruction?
%

% ########################################
% WRITE ANY CODE or COMMENTS HERE to Support your FINAL CONCLUSION
% Is SVD better? If so, why?


% ########################################