%% PCA: Eigen-images

clear variables
clc

%
% You have already done PCA on shape data in a lab to find the elongation and orientation of shapes in 2 spatial dimensions (x- and y-coordinates). Now we'll do PCA on 2+ intensity dimensions.
%% Experiment 1: PCA on three-pixel images
% Consider this simple single-column, 3-pixel image. We'll zoom in to make the pixels look big.

fprintf("A three-pixel image")
img = [.1, .2, .9]'

scale=128
imshow(pixelate(img,scale))

% Plotted as a point in 3D space, where the x value is pixel 1's value, the y value is pixel 2 and the z value is pixel 3, it looks like this. Feel free to rotate the 3D plot.
%

fprintf("The same 3-pixel image. Each dimension shows how light or dark one pixel is")
plot3(img(1), img(2), img(3), 'o')
setLimitsAndLabels(0,1,0,1,0,1,'pix1 intensity', 'pix2 intensity', 'pix3 intensity')

%
% Now, we'll look at a whole collection of 3-pixel images, one per column. In each, the first two pixels are approximately the same intensity and the third pixel is approximately opposite the other two in intensity (gray is its own opposite).

% Set up random number generator with the fixed seed 17 for debugging.
% You can comment out that line later to make it pseudorandom.
rng('default')
rng(17)


noise = 0.1;
unifrnd(-noise, noise, [1,100]);
pixel1 = unifrnd(0, 1, [1,100]);
pixel2 = pixel1 + unifrnd(-noise, noise, [1,100]); % noise added
pixel3 = 1 - pixel1 + unifrnd(-noise, noise, [1,100]); % inverted intensity
imgs = [pixel1; pixel2; pixel3];

% Each column is a 3-pixel image. 
% Show the first few:
fprintf('Five 3-pixel images, one per column')
imshow(pixelate(imgs(:,1:5),scale))

% And plot the whole group of them. 
% In the next step, we'll add the centroid and principal axis. 
fprintf('All the 3-pixel images')
plot3(imgs(1,:), imgs(2,:), imgs(3,:), 'o');
setLimitsAndLabels(0,1,0,1,0,1,'pix1 intensity', 'pix2 intensity', 'pix3 intensity')

%
% It is clear that the values in each dimension are correlated. We perform PCA on this set of points. Complete the findPrincipalComponents() method at the bottom of this file and then call it:

fprintf('Eigvecs, eigvals, and means from findPrincipalComponents')
[V, D, means] = findPrincipalComponents(imgs);

fprintf('Eigenvectors, one per column corresponding to the eigenvalues')
fprintf('The first one is the principal axis')
V

fprintf('All the 3-pixel images again, with the centroid and principal axis added')
hold on
plot3(imgs(1,:), imgs(2,:), imgs(3,:), 'go');
setLimitsAndLabels(0,1,0,1,0,1,'pix1 intensity', 'pix2 intensity', 'pix3 intensity')
plot3([means(1),means(1)+V(1,1)], [means(2),means(2)+V(2,1)], [means(3),means(3)+V(3,1)], 'r-')
plot3([means(1),means(1)-V(1,1)], [means(2),means(2)-V(2,1)], [means(3),means(3)-V(3,1)], 'k-')
plot3(means(1), means(2), means(3), 'r.', 'MarkerSize', 24)
hold off;

fprintf("The principal component is in column 1.")
imshow(pixelate(V,scale))

% How much variation is captured in each eigenvector? There is an easy way to tell: the size of the corresponding eignvalue! (Remember that elongation, ignoring the square root, was the ratio of the eigenvalues.) A long thin shape had a greater ratio, and greater elongation. The eigenvalues returned by eigs will be sorted in descending order. If the first one is much greater than the other two, then it captures almost all of the variability (like a long, thin cluster of points). If they are all close to equal, then the variability is the same in each dimension (like a sphere of points).
% Look at the eigenvalues now.

fprintf('Eigenvalues, sorted from largest to smallest')
diag(D)


% You should find that the first eigenvalue is relatively large. Therefore the first eigenvector (the principal component) captures almost all of the variability - the other two are just capturing slight variations. And indeed, the first column should look like expected: the first 2 pixels are similar and the 3rd is opposite. The second and third columns probably don't make much sense and that's fine - we ignore them since they correspond to eigenvectors with tiny eigenvalues.
%
%% Experiment 2: The three pixels in a different order
% If we swap pixels 2 and 3, the dimensions just change places. Confirm this.

fprintf("Experiment 2: swapping pixels 2 and 3 in each image")
noise = 0.1;
unifrnd(-noise, noise, [1,100]);
pixel1 = unifrnd(0, 1, [1,100]);
pixel2 = 1 - pixel1 + unifrnd(-noise, noise, [1,100]); % Now pixel 2 is inverted
pixel3 = pixel1 + unifrnd(-noise, noise, [1,100]);
imgs = [pixel1; pixel2; pixel3];
imshow(pixelate(imgs(:,1:5),scale))

[V,D, means] = findPrincipalComponents(imgs);
plot3(imgs(1,:), imgs(2,:), imgs(3,:), 'go');
hold on
setLimitsAndLabels(0,1,0,1,0,1,'pix1 intensity', 'pix2 intensity', 'pix3 intensity')
plot3([means(1),means(1)+V(1,1)], [means(2),means(2)+V(2,1)], [means(3),means(3)+V(3,1)], 'r-')
plot3([means(1),means(1)-V(1,1)], [means(2),means(2)-V(2,1)], [means(3),means(3)-V(3,1)], 'k-')
plot3(means(1), means(2), means(3), 'r.', 'MarkerSize', 24)
hold off;
V
D
imshow(pixelate(V,scale))

%% Experiment 3: More pixels
% You have seen that PCA captures correlations between pixels in images. If we move beyond 3 pixels, we won't be able to make nice 3D plots of the point clusters. But the principal components (the eigenvectors corresponding to the max eignvalues), displayed as images, still show which pixels are correlated! What if instead of toy images, we used real images where we expect correlation due to the subject matter? Like a video? Or a sequence of photos taken of the same location every 10 minutes for a whole day, or every day for a year? Or a bunch of Rose-Hulman student ID photos?
% Now is your turn to explore with a bit less direction.
% First, download one of the zipped sets of photos here and unzip them in a convenient location. 4usedinclass is probably the cleanest data set.
% https://www.rose-hulman.edu/class/cs/binaries/csse463/lab7/
%

% TODO: Change to your folder 
imagePath = fileparts('/Users/boutell/Documents/Rose/463/463repo/Public/Exercises/pca/');
% Change if you downloaded/unzipped a different folder 
imageDir = [imagePath, '/4usedinclass'] 

scale2 = 8 % Reduce for now to make it run faster
[imgs, numRows, numCols] = readImageMatrixOnePerColumn(imageDir, scale2);
numRows
numCols

% TODO: Now perform PCA on this matrix like above. Store the top 6 principal 
% eigenvectors and eigenvalues in a matrix called V, each in 1 column, 
% as returned by MATLAB's eig/eigs. You should have already written the
% function to do this for experiment 1.


% Now we want to see what the top eigenvector looks like... as an image: an "eigenimage".
%

% TODO 1: Reshape the first eigenvector to the original dimensions of the
% image. Hint: look at MATLAB's reshape() function and the "reshaping hint 
% #1" below.  
% For the reshape dimensions, note that readImageMatrix...() function 
% returned more than images! And they are color images so need 3 channels. 
% If you still can't figure it out after a couple tries, then
% see the hint at the bottom of the doc.


% TODO 2: Display it pixellated with scale2 since it was scaled down earlier. 
% Note that since eigenvectors have length of 1, you may also 
% have to multiply it by a scale factor (maybe 20 or more), 
% possibly even a negative one to reverse the direction, 
% to get something worth seeing. 



%

% TODO: Repeat for the second and third eigenvectors.


%
%% Experiment 3a: Approximating the image using the principal components
% The final thing to note is that we can approximate the original data using the principal components.
% Let's approximate the image in column 1.
% Step 1: Find the coordinates of this image using the new principal axes. In the following code, I'm assuming that the top two, V1 and V2, capture enough of the variability to be a good approximation. Since V1 and V2 are unit vectors, we just need to take the dot products of the image and each of V1 and V2. I do this with matrix multiplication. This gives us two coordinates in a 2D (plane) that approximates the image set.

% First 2 eigenvectors * the first image. Check the dimensions. The answer
% should be 1x2.
projectionCoordinates = V(:,1:2)'*imgs(:,1) 
size(projectionCoordinates)

%
% This plane is a subspace of the image set, but if it captures most of the variability, then the approximations are close to the original points.
% Step 2: Approximate the image by taking the linear combination of the eigenvectors and then adding back the means. If the projection coordinates are denoted as pc, then
% approximation = means +  sum from i = 1 to 2 of pc(i)*V(i)

% TODO: Calculate the approximation using the linear combination and means.  
% Then reshape the result back to the image shape. 
% Then display it pixellated. 
% You should get something that looks like the first image in the set.


%
% Reshaping hint #1:

% How the colon operator works on 3D matrices:
a(:,:,1) = [1 2; 3 4];
a(:,:,2) = [5 6; 7 8];
a
size(a)
b = a(:)

% And how reshape undoes them. 
c = reshape(b, 2,2,2)

%
% Complete this method, to be used in all the experiments. It's mostly done for you.

function [V, D, means]= findPrincipalComponents(data)
    % Input: data is a dataset, one example per column. Each column 
    % could be a whole image, for example.
    % size(data, 1) is the dimensionality of each example.
    % size(data, 2) is the number of examples in the data set. 
    % Outputs:
    % V: the first principal components of the dataset, found by taking the 
    % eigenvectors of the data's covariance matrix. Return at most 6 of
    % them.
    % D: the corresponding eigenvalues of the covariance matrix, stored
    % along the diagonal of the D matrix.
    % means: the centroid of the data set.


    % Shift the data cloud so that its centroid is at the origin.
    means = mean(data, 2); % mean of each row
    normed = data - means;
    
    % You are given the means, and the other 2 return values seem very
    % specific! But don't worry. Read the MATLAB docs for the EIGS function
    % (which is more efficient than EIG with large matrices when you only 
    % want the most important eigenvectors) and find the version that 
    % returns V and D. Then you'll see you only have a little more code
    % to write in this function! It can be a one-liner.
    % TODO
    
    
end 

% Given helper function(s):
% Make sure you understand how we are storing the images in the matrix.

function [imgs, numRows, numCols] = readImageMatrixOnePerColumn(imagedir, scale)
    fileList = dir(fullfile(imagedir,'*.jpg'));
    
    % Get the number of images and the number of pixels in each image 
    % so that we can pre-allocate a matrix to hold all the images, 
    % one per column.
    % Each image MUST have the same dimensions!
    nImages = size(fileList, 1) - 2; % -2 to remove the . and .. folders.

    % Get the dimension of 1 image
    img = imread([imagedir  '/'  fileList(3).name]);
    img = imresize(img, 0.125, 'bicubic');
    numRows = size(img, 1);
    numCols = size(img, 2);
    nPixels = size(img(:),1);
    imgs = zeros(nPixels, nImages);
    
    % files 1 and 2 are . (current dir) and .. (parent dir), respectively, 
    % so we start with 3.
    % We reduce each image by a scale factor so it runs faster and 
    % uses less memory.
    for i = 3:size(fileList)
        img = imread([imagedir  '/'  fileList(i).name]);
        img = imresize(img, 1/scale, 'bicubic');
        % The most important line. Understand it?
        imgs(:, i-2) = img(:);
    end
    % Convert to [0,1] double. 
    imgs = double(imgs);
    imgs = imgs / 255; 
end 


function out = pixelate(img, scale)
    [nr,nc,nchannels] = size(img);
    out = zeros(nr*scale, nc*scale, nchannels);
    for r = 1:nr
        for c = 1:nc
            pixel = img(r,c,:);
            replicatedPixel = repmat(pixel, [scale, scale, 1]);
            % size(replicatedPixel);
            out((r-1)*scale+1:r*scale, (c-1)*scale+1:c*scale, :) = replicatedPixel;
        end
    end
end

function setLimitsAndLabels(xmin, xmax, ymin, ymax, zmin, zmax, xlab, ylab, zlab)
    xlim([xmin xmax])
    ylim([ymin ymax])
    zlim([zmin zmax])
    xlabel(xlab) 
    ylabel(ylab) 
    zlabel(zlab) 
end


% Still need a hint? Here is a key line of reshaping code used in a couple places.

% reshapedImage = reshape(columnImage, numRows, numCols, numChannels);
% numChannels is 3 since RGB.