% First Implementation Attempting to Define Nucleus Space
% May 24, 2007
% (c) 2007 Christopher Mitchell
% All rights reserved

disp('Stage 1 (reading training images)...');
% First we need to load a bunch of images from file.
% Note: this program calculates dynamic array sizes, but it assumes
%       that all images have the same (square) dimensions.  Make it so!
RGB1 = imread('./v1/subject01.normal.png');
RGB2 = imread('./v1/subject02.normal.png');
RGB3 = imread('./v1/subject03.normal.png');
RGB4 = imread('./v1/subject04.normal.png');
RGB5 = imread('./v1/subject05.normal.png');
RGB6 = imread('./v1/subject06.normal.png');
RGB7 = imread('./v1/subject07.normal.png');
RGB8 = imread('./v1/subject08.normal.png');

disp('Stage 2 (generating matrix from training set)...');
% Next we need to create a big matrix containing all of the data we just
% loaded as 8 matrices
[dimY,dimX] = size(RGB1);
RawData = zeros(dimY*dimX,dimY*dimX);
for i = 1:dimY
    for j = 1:dimX
        RawData(1,((i-1)*dimX)+j) = RGB1(i,j);
        RawData(2,((i-1)*dimX)+j) = RGB2(i,j);
        RawData(3,((i-1)*dimX)+j) = RGB3(i,j);
        RawData(4,((i-1)*dimX)+j) = RGB4(i,j);
        RawData(5,((i-1)*dimX)+j) = RGB5(i,j);
        RawData(6,((i-1)*dimX)+j) = RGB6(i,j);
        RawData(7,((i-1)*dimX)+j) = RGB7(i,j);
        RawData(8,((i-1)*dimX)+j) = RGB8(i,j);
    end
end

disp('Stage 3 (calculating eigenCells)...');
% Step 3: Compute and sort eigenvectors and eigenvalues
[eigenVectors,eigenValues] = eig(RawData);
eigenValues = diag(eigenValues);
eigenVectors(:,dimX*dimY+1) = eigenValues;
eigenVectors = sortrows(eigenVectors,dimX*dimY+1);
eigenVectors = eigenVectors(:,1:dimX*dimY);
eigenVectors = eigenVectors(1:8,:);

disp('Stage 4 (process test image)...');
% first, input a file and pack into a vector
file =  input('Enter a filename to test:','s');
RGB9 =  imread(file);
RGB9v = zeros(1,dimY*dimX);
for i = 1:dimY
    for j = 1:dimX
        RGB9v(1,((i-1)*dimX)+j) = RGB9(i,j);
    end
end

% now, transform from image space to face space
cellMatch = zeros(1,8);
cellMatch(1,1) = inner(RGB9v,eigenVectors(1,:));
cellMatch(1,2) = inner(RGB9v,eigenVectors(2,:));
cellMatch(1,3) = inner(RGB9v,eigenVectors(3,:));
cellMatch(1,4) = inner(RGB9v,eigenVectors(4,:));
cellMatch(1,5) = inner(RGB9v,eigenVectors(5,:));
cellMatch(1,6) = inner(RGB9v,eigenVectors(6,:));
cellMatch(1,7) = inner(RGB9v,eigenVectors(7,:));
cellMatch(1,8) = inner(RGB9v,eigenVectors(8,:));

% next, convert back into image space
RGB9v2 = zeros(1,dimY*dimX);
RGB9v2 = RGB9v2 + cellMatch(1,1)*eigenVectors(1,:);
RGB9v2 = RGB9v2 + cellMatch(1,2)*eigenVectors(2,:);
RGB9v2 = RGB9v2 + cellMatch(1,3)*eigenVectors(3,:);
RGB9v2 = RGB9v2 + cellMatch(1,4)*eigenVectors(4,:);
RGB9v2 = RGB9v2 + cellMatch(1,5)*eigenVectors(5,:);
RGB9v2 = RGB9v2 + cellMatch(1,6)*eigenVectors(6,:);
RGB9v2 = RGB9v2 + cellMatch(1,7)*eigenVectors(7,:);
RGB9v2 = RGB9v2 + cellMatch(1,8)*eigenVectors(8,:);

%finally, find the mean-squared error between original and this image
MSE = sum((RGB9v2-RGB9v).^2);
disp('Mean squared error:');
disp(MSE);
RGB92 = zeros(dimY,dimX);
for i = 1:dimY
    for j = 1:dimX
        RGB92(i,j) = RGB9v2(1,j+dimX*(i-1));
    end
end
figure
imshow(RGB9);
figure
imshow((RGB92-min(min(RGB92)))/(max(max(RGB92))-min(min(RGB92))));
disp('Program termination (normal): completed.');