% Third Implementation (Search Subprogram)
% May 29-31, 2007
% (c) 2007 Christopher Mitchell
% All rights reserved

% constants
threshold = 4e7;

% First we need to load a bunch of images from file.
% Load them from the preexisting training file
load eigenCells eigenCells
eigenSize = size(eigenCells);
eigenN = eigenSize(1);
eigenX = input('Width of training images (px): ');
eigenY = input('Height of training images (px): ');
if eigenX*eigenY ~= eigenSize(2)
    disp('Dimensional mismatch in training set.');
    return  %This would be a dim mismatch
end

% Set up matrices beforehand for speed
s = input('Dimensions match.  Now we need the test image: ','s');
testImage = imread(s);
[IdimY,IdimX] = size(testImage);
matches = zeros(IdimY-eigenY+1,IdimX-eigenX+1);
subimage = zeros(eigenY,eigenX);
subimagevector = zeros(1,eigenX*eigenY);
cellMatch = zeros(eigenN);

% Now do the tests for each possible subimage.  This should be a LONG
% process since we're operating on a 250,000-element matrix. :P
tic;
for X = 1:size(matches,2)
    for Y = 1:size(matches,1)
        subimage = testImage(Y:Y+eigenY-1,X:X+eigenX-1);
        for i = 1:eigenY
            for j = 1:eigenX
                subimagevector(1,((i-1)*eigenX)+j) = subimage(i,j);
            end
        end
        %now do the facespace transformations
        
        %% Commented-out percentage printing code (needs buffer flushing) %%
        %if floor((X*size(matches,1)+Y)/1000) == (X*size(matches,1)+Y)/1000
        %    s = sprintf('%d%% complete',floor(100*((X*size(matches,1)+Y)/(size(matches,1)*size(matches,2)))));
        %    disp(s);
        %end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        %To face space
        for i=1:eigenN
            cellMatch(i) = inner(subimagevector,eigenCells(i,:));
        end
        
        %Back to image space
        subimagevector2 = zeros(1,eigenX*eigenY);
        for i=1:eigenN
            subimagevector2 = subimagevector2 + cellMatch(i)*eigenCells(i,:);
        end
        
        %Find the MSE
        matches(Y,X) = sum(((subimagevector2-min(min(subimagevector2)))/(max(max(subimagevector2))-min(min(subimagevector2)))-subimagevector).^2);
    end
end
toc;

matchesBox = ones(IdimY-eigenY+1,IdimX-eigenX+1);
maxval = max(max(matches));
while min(min(matches)) < threshold
    [Y,I] = min(matches);
    [X,J] = min(Y);
    coordX = J;
    coordY = I(J);
    tic;
    s = sprintf('(%d,%d)',coordX,coordY);
    disp(s);
    for X = coordX-floor(eigenX/2):coordX+floor(eigenX/2)
        for Y = coordY-floor(eigenY/2):coordY+floor(eigenY/2)
            if X>0 && Y>0 && X<=size(matches,2) && Y<=size(matches,1)
               matches(Y,X) = maxval;
               matchesBox(Y,X) = 0;
               if X == coordX-floor(eigenX/2) || X == coordX+floor(eigenX/2) || Y == coordY-floor(eigenY/2) || Y == coordY+floor(eigenY/2)
                   matchesBox(Y,X) = 1;
               end
            end
        end
    end
    toc;
end

matches = matches-min(min(matches));
matches = matches./(max(max(matches)) - min(min(matches)));
imshow(matches);

%All done here!