function destImg = findEye(sourceImg)
% destImg = findeye(sourceImg).  This function reads in an image and marks
% the eyes it locates within that image.
% Author : Sasakthi Abeysinghe
% Date   : 01/20/2006

global imageDirectory;

destImg = sourceImg;
sourceImg = double(rgb2gray(sourceImg)) / 256.0;

height = size(sourceImg, 1);
width = size(sourceImg, 2);
% Load the template Eye.  This should work better if the template eye was
% an averaged image of a lot of eyes, with an additional matrix with the
% standard deviation which should be used to weigh the pixel locations.
template = double(rgb2gray(imread(strcat(imageDirectory, 'Eye.jpg')))) / 256.0;

% The number of iterations are calculated assuming that the size of eye
% will not exceed 25% of the size of the image.
iterations = ceil(width / 40);

% A matrix storing the ssd scores between the image and the template.
% Initialized at a very high value, but the actuall values will be between
% 0 and 1.
ssd = ones(height, width, iterations) * 1000;

% Looping to try different sizes of eyes.
for i = 1:iterations,
    w = i * 10;
    eyeMask = imresize(template, w / size(template, 2));
    maskHeight = size(eyeMask, 1);
    maskWidth = size(eyeMask, 2);
    normalizingFactor = maskHeight * maskWidth;
    
    for x = 1:width-maskWidth,
        for y = 1:height - maskHeight,
            imageMask = sourceImg(y : y + maskHeight - 1, x : x + maskWidth - 1);
            delta = eyeMask - imageMask;           
            ssd(y + ceil(maskHeight / 2), ceil(x + maskWidth / 2), i) = sum(sum(delta .* delta)) / normalizingFactor;            
        end
    end
end

threshold = min(min(min(ssd))) + 0.015;
for r = 1:height,
    for c = 1:width,
        for z = 1:iterations,
            if (ssd(r, c, z) < threshold)
                r1 = max(floor(r-z*2.5), 1);
                r2 = min(floor(r+z*2.5), height);
                c1 = max(c-z*5, 1);
                c2 = min(c+z*5, width);
                % Marking a box around the identified eye location.
                destImg(r1, c1:c2, 1:3) = 255;
                destImg(r2, c1:c2, 1:3) = 255;
                destImg(r1:r2, c1, 1:3) = 255;
                destImg(r1:r2, c2, 1:3) = 255;       
                % Clearing other potential eye locations in the local
                % vicinity
                ssd(r1:r2, c1:c2, 1:iterations) = threshold + 1;
            end
        end
    end
end
image(destImg);