% Eric La Force
% COMP 776 - Assignment 01
% Blob Detector
function detectblobs(filename, s, scaleFactor, thres, num, method)
tic
%% INITIALIZATION
image = im2double(rgb2gray(imread(filename)));
numberOfIterations = num;
sigma = s;
k = scaleFactor;
th = thres;
numberOfCircles = 0;
cx = [];
cy = [];
cr = [];
sigmaValues = [];
if(strcmp(method,'brute') == 1)
    filter = fspecial('log',30,sigma);
else
    filter = sigma^2*fspecial('log',30,sigma);
end

[heightOfImage, widthOfImage ] = size(image);
scaleSpace = zeros(heightOfImage,widthOfImage,numberOfIterations); % [h,w] - dimensions of image, n - number of levels in scale space
maxSpace = zeros(heightOfImage,widthOfImage,numberOfIterations);

%% Create Scale Space
temp = image;
for i=1:numberOfIterations 
    sigmaValues(i) = sigma;
    if(strcmp(method,'brute') == 1)
        normalizedFilter = imfilter(image,sigma^2*filter,'replicate');
    else
        temp = imresize(temp,1/k,'bicubic');
        normalizedFilter = imfilter(temp,filter,'replicate');
    end
    
    squaredNormalizedFilter = normalizedFilter.^2;
    
    if(strcmp(method,'brute') == 0)
        squaredNormalizedFilter = imresize(squaredNormalizedFilter,[heightOfImage widthOfImage],'bicubic');
    end
    
    scaleSpace(:,:,i) = squaredNormalizedFilter(:,:);
    f = ordfilt2(squaredNormalizedFilter,25,true(5));
    maxSpace(:,:,i) = f(:,:);
    sigma = sigma*k;
    
    if(strcmp(method,'brute') == 1)
        filter = fspecial('log',30,sigma);
    end
end

%% Preform non maximum suppression for the scale space
 for i=1:numberOfIterations
    for h=floor(sigma):heightOfImage - floor(sigma)
       for w=floor(sigma):widthOfImage - floor(sigma)
           pix = scaleSpace(h,w,i);
           if( (pix > th) && ...
               (pix > scaleSpace(h-1,w-1,i)) && ...
               (pix > scaleSpace(h,w-1,i)) && ...
               (pix > scaleSpace(h+1,w-1,i)) && ...
               (pix > scaleSpace(h-1,w,i)) && ...
               (pix > scaleSpace(h+1,w,i)) && ...
               (pix > scaleSpace(h-1,w+1,i)) && ...
               (pix > scaleSpace(h,w+1,i)) && ...
               (pix > scaleSpace(h+1,w+1,i)) && ...
               (((i+1) < 11) && pix > maxSpace(h,w,i+1)) && ...
               (((i-1) > 0) && pix > maxSpace(h,w,i-1)))
                  numberOfCircles = numberOfCircles + 1;
                  cx(1,numberOfCircles) = w;
                  cy(1,numberOfCircles) = h;
                  cr(1,numberOfCircles) = sigmaValues(i)*sqrt(2);
                  
           end
       end
    end
end


%% 4. Display resulting circles at their characteristic scales. 
show_all_circles(image,cx',cy',cr','g');
toc