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

% 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): ');
threshold = input('Threshold: ');
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');
blockX = input('Density block width: ');
blockY = input('Density block height: ');

testImage = imread(s);
[IdimY,IdimX] = size(testImage);
density = zeros(ceil(IdimY/blockY),ceil(IdimX/blockX));
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);
%        if max(max(subimage)) ~= min(min(subimage))
%            subimage = subimage./(max(max(subimage))-min(min(subimage)));
%        end
%        subimage = subimage - min(min(subimage));
        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
        if (max(max(subimagevector2))~=min(min(subimagevector2)) && (sum(((subimagevector2-min(min(subimagevector2)))))~=0))
            matches(Y,X) = sum(((subimagevector2-min(min(subimagevector2)))/(max(max(subimagevector2))-min(min(subimagevector2)))-subimagevector).^2);
        else
            matches(Y,X) = threshold.^2;
        end
    end
end
toc;

matchesBox = ones(IdimY-eigenY+1,IdimX-eigenX+1);
maxval = max(max(matches));

RGBimg = zeros(size(testImage,1),size(testImage,2),3);
RGBimg(:,:,1) = testImage;
RGBimg(:,:,2) = testImage;
RGBimg(:,:,3) = testImage;

while min(min(matches)) < threshold
    [Y,I] = min(matches);
    [X,J] = min(Y);
    coordX = J;
    coordY = I(J);
    s = sprintf('(%d,%d) @ MSE threshold of %d',coordX,coordY,X);
    disp(s);
    density(ceil(coordY/blockY),ceil(coordX/blockX)) = density(ceil(coordY/blockY),ceil(coordX/blockX))+1;
    for X = coordX:coordX+eigenX-1
        for Y = coordY:coordY+eigenY-1
            if (X>0 && Y>0 && X<=size(RGBimg,2) && Y<=size(RGBimg,1)) && (X == coordX || X == coordX+eigenX-1 || Y == coordY || Y == coordY+eigenY-1)
               RGBimg(Y,X,2) = 255;
               RGBimg(Y,X,1) = 0;
               RGBimg(Y,X,3) = 0;
            end
        end
    end
    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) = threshold.^2;
            end
        end
    end
end

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

[densityimy,densityimx] = size(density);
r = zeros(densityimy,densityimx);
g = r;
b = r;
density = density - min(min(density));
density = density/max(max(density));
r = 2*density;
g = 2*(1-density);
figure;

% now it's time to expand it up to full size! %
RGBimgD = zeros(size(testImage,1),size(testImage,2),3);
for x=1:size(testImage,2)
    for y=1:size(testImage,1)
        RGBimgD(y,x,1) = r(ceil(y/blockY),ceil(x/blockX));
        RGBimgD(y,x,2) = g(ceil(y/blockY),ceil(x/blockX));
        RGBimgD(y,x,3) = b(ceil(y/blockY),ceil(x/blockX));
    end
end
imshow(RGBimgD);

%All done here!