% Third Implementation (Training Set Subprogram)
% May 26-31, 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!
count = input('How many training images? ');
dimX = input('Image width (px)? ');
dimY = input('Image height (px)? ');

grays = zeros(dimX,dimY);
RawData = zeros(count,dimY*dimX);

disp('Stage 2 (generating matrix from training set)...');
for k = 1:count
    s = input('Filename:','s');
    grays = imread(s);
    grays = grays./(max(max(grays))-min(min(grays)));
    grays = grays - min(min(grays));
    for i = 1:dimY
        for j = 1:dimX
            RawData(k,((i-1)*dimX)+j) = grays(i,j);
        end
    end
end

disp('Stage 3 (calculating outer product)...');
% Step 3: Calculate the count*count matrix
outerMat = zeros(count,count);

for i=1:count
    for j=1:count
        outerMat(i,j) = inner(RawData(i,:,:),RawData(j,:,:));
    end
end
[eigenVectors,eigenValues] = eig(outerMat);
eigenValues = diag(eigenValues);

disp('Stage 4 (generating eigenCells)...');
% Step 4: Use the eigenvectors of the outer matrix to generated eigenCells
eigenCells = zeros(count,dimX*dimY);
eigenCellImages = zeros(count,dimX,dimY);
for i=1:count
    for j=1:count
        eigenCells(i,:) = eigenCells(i,:) + eigenVectors(i,j).*RawData(j,:);
    end
end

save eigenCells eigenCells
disp('Program termination (normal): completed.');