function scale_space = blob(imname, start_perc, end_perc, steps)
%% imname     - the filename to read
%% start_perc - the starting percent of the kernel to the image size
%% end_perc   - the ending percent of the kernel to the image size
%% steps      - the number of steps between the start and end percent
%%
%% A sample call:
%% >> blob('butterfly.jpg', 1.5,20,20);
%%    reads 'butterfly.jpg', samples with kernel size ranging from 1.5% 
%%    to 20% of image size, in 20 steps (plus 2 extra samplings)
    img = imread(imname);

    % convert to gray.
    img_float = double(rgb2gray(img));
    image_size = size(img_float);
    scale_space = zeros(image_size(1), image_size(2), steps);
    max_results = zeros(image_size(1), image_size(2), steps);

    % make the kernel
    sigma = 3;
    kern_size = max(1,fix(6*sigma));    % as harris set it.
    kern = fspecial( 'log', kern_size, sigma );
    norm_kern = sigma.^2 * kern; % normalized it.
    %figure(1);
    %imagesc(norm_kern);
    
    % the uppermost and lowermost scales should not be tested for maximums
    % so add in 2 more sample steps
    steps = steps + 2;


    % apply the filter for each level
    for i=1 : steps
        %resize, filter, scale back
        new_image_size = kern_size / percent_at_level(i, start_perc, end_perc, steps);
        resized_img = imresize( img_float, new_image_size / image_size(2), 'bilinear');
        filtered_img = imfilter( resized_img, norm_kern, 'replicate');
        filtered_img = filtered_img .^2;
        filtered_img = imresize( filtered_img, [image_size(1) image_size(2)],  'bilinear');
        
        %dump info
        percent = [num2str(  percent_at_level(i, start_perc, end_perc, steps)*100 ) '%'];
        start_size = [ num2str(image_size(1)) 'x' num2str(image_size(2))];
        end_size = [ num2str( fix(new_image_size / image_size(2) * image_size(1))) 'x' num2str( fix(new_image_size / image_size(2) * image_size(2)))];
        [' '  percent ' = ' start_size ' -> ' end_size]

        scale_space(:,:, i) = filtered_img;
    end

    % 2d nonmax suppression, taken from harris code
    possible_max_list = [];
    for i=1 : steps
        % match the dilated image and are also greater than the threshold.
        current_image = scale_space(:,:,i);
        %threshold = median(current_image(:))
        threshold = 1000;
        %mask_size = fix(image_size(2) * (start_perc + (i-1) * step_size) / 1000);
        mask_size = 5;
        mx = ordfilt2(scale_space(:,:,i), mask_size^2, ones(mask_size)); % Grey-scale dilate.
        max_results(:,:,i) = (scale_space(:,:,i)==mx)&(scale_space(:,:,i)>threshold);       % Find maxima.

        [r,c] = find(max_results(:,:,i));                  % Find row,col coords.
        possible_max_list = [possible_max_list, {[r,c]}];
    end
    
    
    %find max in scale space - 3d nonmax suppression
    max_x = [];
    max_y = [];
    max_r = [];
    for level=2 : steps - 1
        current_max_level = possible_max_list(level);
        current_max_cell = current_max_level{1};
        num_found = size( current_max_cell );
        
        current_radius = fix(image_size(2) * percent_at_level(level, start_perc, end_perc, steps)  * 0.4);
        %filter_size = (5-1) / 2;
        
        % instead of a static filter, scale with relative size of kernel
        filter_size = fix(image_size(2) * percent_at_level(level, start_perc, end_perc, steps)  * 0.5);
        
        ['level: ' num2str(level) ' radius: ' num2str(current_radius) ' filter: ' num2str(filter_size)]
        
        for i=1 : num_found(1)
            current_point = current_max_cell(i,:);
            x = current_point(1);
            y = current_point(2);
            
            % drop the max positions that are too close to the edge of the image
            if(x <= filter_size || x > image_size(1)-filter_size)
                continue;
            end
            if(y <= filter_size || y > image_size(2)-filter_size)
                continue;
            end

            local_max = scale_space(x,y,level);

            upper_vals = [ scale_space(x-filter_size:x+filter_size, y-filter_size:y+filter_size, level+1) ];
            upper_max = max(upper_vals(:));
            
            lower_vals = [ scale_space(x-filter_size:x+filter_size, y-filter_size:y+filter_size, level-1) ];
            lower_max = max(lower_vals(:));
            
            if( (local_max > upper_max) && (local_max > lower_max))
                max_x = [max_x, [x]];
                max_y = [max_y, [y]];
                max_r = [max_r, [current_radius]];
            end
        end
    end
    
    %{
    % display the image
    for level=1 : steps    
        figure(level);
        %imshow( scale_space(:,:, level) );
        subplot(1,2,1);
        imagesc( scale_space(:,:, level) );
        subplot(1,2,2);
        imagesc( max_results(:,:, level) );
    end
    %}
    
    figure(steps+1);
    %show_all_circles(img, max_y', max_x', max_r');%, color, ln_wid);
%    show_circles(img, max_y', max_x', max_r', [0 0 1], 1.5,
    [num2str(start_perc) '% -> ' num2str(end_perc) '% , ' num2str(steps) ' steps']);

end


function perc = percent_at_level(i, start_perc, end_perc, steps)
    step_size = (end_perc - start_perc) / (steps-1);
    % linear steps don't work well
    % instead, sample more often when kernel is small relative to image size
    perc = (((i-1) * step_size) / (end_perc - start_perc))^2 * (end_perc - start_perc) + start_perc;
    %perc = (start_perc + (i-1) * step_size);
    perc = perc / 100;
end

function show_circles(I, cx, cy, rad, color, ln_wid, text)
%% I: image on top of which you want to display the circles
%% cx, cy: column vectors with x and y coordinates of circle centers
%% rad: column vector with radii of circles. 
%% The sizes of cx, cy, and rad must all be the same
%% color: optional parameter specifying the color of the circles
%%        to be displayed (red by default)
%% ln_wid: line width of circles (optional, 1.5 by default

if nargin < 5
    color = 'r';
end

if nargin < 6
   ln_wid = 1.5;
end

imshow(I); hold on;

theta = 0:0.1:(2*pi+0.1);
cx1 = cx(:,ones(size(theta)));
cy1 = cy(:,ones(size(theta)));
rad1 = rad(:,ones(size(theta)));
theta = theta(ones(size(cx1,1),1),:);
X = cx1+cos(theta).*rad1;
Y = cy1+sin(theta).*rad1;
line(X', Y', 'Color', color, 'LineWidth', ln_wid);

%title(sprintf('%d circles', size(cx,1)));
title(text);
end