%read the image

image_name='onecircle.png';
%image_name='scalespaceTest.png';
%image_name='sunflower3.png';
I=imread(image_name);

I = double(I);
im_size=size(I);
%imshow(I);

%create the gaussian scale space and do the laplacian filter
t_array=1:0.1:3;
sigma_array = exp(t_array);
blurred_image = cell(length(t_array),1); %creat a cell, which will contain all the bluured_images
lap_image = cell(length(t_array),1); %cell for the image after laplacian filter
%lap_filter = fspecial('laplacian',0.001); % create a laplacian filter
lap_filter = [0, 1, 0;1, -4, 1; 0, 1, 0];
coef_gau_filter=1.8; %control the size of the gaussian filter 
coef_peak_mask=9;  %control the size of the neighbor



for k = 1:length(t_array)
    gau_filter = fspecial('gaussian', round(sigma_array(k)^coef_gau_filter),sigma_array(k));
    blurred_image{k,1}=imfilter(I, gau_filter,'replicate');
    
    %filter each scale by a discrete Laplacian kernel
    temp_lap_filter = lap_filter;
    lap_image{k,1} = imfilter(blurred_image{k,1},temp_lap_filter, 'replicate') * sigma_array(k)^2;
    
    %save the image
    if k <10
        imwrite(blurred_image{k,1},strcat('.\resultImages\blurred\blurred_image0', num2str(k),'.png'),'png');
        imwrite(lap_image{k,1},strcat('.\resultImages\lap\lap_image0', num2str(k),'.png'),'png');
    else
        imwrite(blurred_image{k,1},strcat('.\resultImages\blurred\blurred_image', num2str(k),'.png'),'png');
        imwrite(lap_image{k,1},strcat('.\resultImages\lap\lap_image', num2str(k),'.png'),'png');
    end
    %k
    %lap_image{k,1}(45,44)
    %figure;
    %imshow(lap_image{k,1});
    
    
end

%find the center points:

%peak_list_lower=[]; % strore the peak for a lower scale
%peak_list_upper=[]; % strore the peak for a uper scale
peak_list=[];
total_list=[];
min_dis=0;
max_intensity=0;
index_del=[];
index_add=[];


for k = 1:length(t_array)
    
    if k ==1
        for m = 1:im_size(1)
            for n = 1: im_size(2)
                if ispeak(lap_image{k,1}, m, n, coef_peak_mask) == 1
                    total_list= [total_list; m,n,k,lap_image{k,1}(m,n)];
                end
            end
        end
    else
        %computhe the upper
        for m = 1:im_size(1)
            for n = 1: im_size(2)
                if ispeak(lap_image{k,1}, m, n, coef_peak_mask) == 1
                    peak_list= [peak_list; m,n,k,lap_image{k,1}(m,n)];
                end
            end
        end
        
        
        
        
        %test if the lower peak is inside of the circles centering at the
        %upper peak
        peak_size=size(peak_list);
        total_size=size(total_list);
        
        if total_size(1) == 0
            total_list=peak_list;
            peak_list=[];
        else
        
        
        
            for p = 1:peak_size(1)
                for q = 1: total_size(1)
                    distance = sqrt(dot(total_list(q,1:2) - peak_list(p,1:2), total_list(q,1:2) - peak_list(p,1:2)));

                    if distance < min_dis || q==1 
                        min_dis = distance;
                        min_index = q;
                    end
                    
                    
                    if distance < sigma_array(k)*1.5
                       if  lap_image{k,1}(peak_list(p,1), peak_list(p,2)) > lap_image{total_list(q,3),1}(total_list(q,1), total_list(q,2))
                          %index_del=add_to_list(index_del,q);
                          index_del=[index_del,q];
                          %index_add=[index_add,p];
                          %total_list = [total_list;peak_list(p,:)];
                       end
                    %else
                    %    index_add=[index_add,p];
                        %total_list = [total_list;peak_list(p,:)];
                    end

                end
                
                
                if min_dis < sigma_array(k)*1.5
                    if  lap_image{k,1}(peak_list(p,1), peak_list(p,2)) > lap_image{total_list(min_index,3),1}(total_list(min_index,1), total_list(min_index,2))
                          %index_del=add_to_list(index_del,q);
                          index_del=[index_del,q];
                          index_add=[index_add,p];
                    end
                else
                    index_add=[index_add,p];   
                end
                
                
               
                
            end

            total_list=delete_from_list(total_list,index_del);
            total_list=add_from_list(total_list, peak_list,index_add);
            index_del=[];
            index_add=[];
            peak_list=[];
            
            %total_list
            %k
        end
        
        
    end
    
    
    
end
total_list
size(total_list)




%draw the circle
drawCircle2(image_name, total_list(:,1),total_list(:,2),sigma_array(total_list(:,3))*1.5);



