% The Cellular Density Project (CDP)
% Final Evaluation Program
% July 13, 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

%default values here
default_image_threshold = 2e7;
%default values end here

%Find out how many eigenCells we have stored from the training program
eigenSize = size(eigenCells);
eigenN = eigenSize(1);

%Verify eigenCell dimensions
eigenX = input('Width of training images (px): ');
eigenY = input('Height of training images (px): ');
s = sprintf('Threshold (default=%d): ',default_image_threshold);
threshold = input(s);
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: ');
disp(' ');

% Read in the image to be processed
testImage = imread(s);
[IdimY,IdimX] = size(testImage);

% Create the density array for use later
density = zeros(ceil(IdimY/blockY),ceil(IdimX/blockX));
% Create the matches (MSEs values) for use later
matches = zeros(IdimY-eigenY+1,IdimX-eigenX+1);
% All the subimages are the same size, so we can pre-created that too
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)
        % Create this subimage
        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
        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
            % Error-catching: set to a high MSE if we would encounter a div
            % by zero situation
            matches(Y,X) = threshold.^2;
        end
    end
end
toc;

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

% Set up the RGB image (grayscale with green boxes)
RGBimg = zeros(size(testImage,1),size(testImage,2),3);
RGBimg(:,:,1) = testImage;
RGBimg(:,:,2) = testImage;
RGBimg(:,:,3) = testImage;

% This nested bunch of stuff finds all matching cells below the
% threshold for cell recognition
nummatches = 0;
while min(min(matches)) < threshold
    [Y,I] = min(matches);
    [X,J] = min(Y);
    coordX = J;
    coordY = I(J);
    nummatches = nummatches+1;
    % Extra debug code commented out
    % 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;
    % The following code erases the area of the cell from the MSE array and
    % adds the green box to the RGB image
    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

% Some output so the user knows how far we've gotten
s = sprintf('%d cells were matched in the image.',nummatches);
disp(s);
disp('Generating hitmap...');

% Show the green-boxed match image
matches = matches-min(min(matches));
matches = matches./(max(max(matches)) - min(min(matches)));
imshow(RGBimg./255);

disp('Generating density heatmap...');
[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
% Show the density heatmap image
imshow(RGBimgD);

% Now do all the statistical analysis stuff
disp('Evaluating cancer presence in normalized density map');
dstd = std(density(:));
dmean = mean(density(:));
dmin = min(min(density));
dmax = max(max(density));
s = sprintf('Mean of %1.4d and standard deviation of %1.4d found.',dmean,dstd);
disp(s);
dminpct = 100*((dmax-dmean)/dstd);
dmaxpct = 100*((dmean-dmin)/dstd);
s = sprintf('Min is %4.1d%% std dev below mean',dminpct);
disp(s);
s = sprintf('Max is %4.1d%% std dev above mean',dmaxpct);
disp(s);

% Find out how many are inside and out, and the averages in terms of the
% std dev
inside = 0;
insidediff = 0;
outside = 0;
outsidediff = 0;
for x=1:size(density,2)
    for y=1:size(density,1)
        if ((density(y,x) < dmean+dstd) && (density(y,x) > dmean-dstd))
            inside = inside+1;
            insidediff = insidediff + (abs(density(y,x)-dmean)/dstd);
        else
            outside = outside+1;
            outsidediff = outsidediff + (abs(density(y,x)-dmean)/dstd);
        end
    end
end

% Final analysis output
s = sprintf('%d regions inside, %d regions outside',inside,outside);
disp(s);
s = sprintf('Inside stdev average at %4.1d%%, outside stddev average at %4.1d%%',100*(insidediff/inside),100*(outsidediff/outside));
disp(s);
disp(' ');
uncertainty = 100*min([inside/outside outside/inside]);
if (inside > outside)
    s = sprintf('*Based on deviation:* This image does not contain cancerous region(s) at uncertainty of %d%%.',uncertainty);
else
    s = sprintf('*Based on deviation:* This image contains cancerous region(s) at uncertainty of %d%%.',uncertainty);
end
disp(s);
if (0.4 < dmean && dmean < 0.6)
    s = sprintf('*Based on mean:* This image does not contain cancerous regions at uncertainty of %4.1d%%.',200*abs(dmean-0.5));
else
    s = sprintf('*Based on mean:* This image contains cancerous region(s) at uncertainty of %4.1d%%.',200*(0.5-abs(dmean-0.5)));
end
disp(s);
disp(' ');
disp('Program termination: program complete.');
disp(' ');

%All done here!