%------------------------------------------------------------------
%
% edge_script.m
%
% Script that demonstrates the process of detecting edges in an image
% This is essentially the canny edge detector
%
%-------------------------------------------------------------------



%%-- Image name:

name_image_list = {'fireplace.bmp', 'test_ima.bmp'}

ima_number = input('Which image do you want to process (1=fireplace.bmp, 2=test_ima.bmp) ? ');

name_image = name_image_list{ima_number};


%%-- Load image:

ima = imread(name_image);

%%-- Display the image:

figure(1);
image(ima);
title('Input image');
axis image
drawnow;

% Grayscale convertion:
ima2 = rgb2gray(ima);

figure(2);
image(ima2);
colormap(gray(256));
title('Input image - (converted in B&W)');
axis image
drawnow;

% Convert i double (from uint8)
ima3 = double(ima2);

[m,n] = size(ima3);


%%%-- Parameters of edge detection:
hithresh = 0.7;
fitthresh = 0.1;
lothresh = hithresh/2;
scale = 1;

%% construct blurring kernel k (it is a 1D Gaussian)
halfklen = 3 * scale;
ksupport = -halfklen:1:halfklen;

k = exp( - (ksupport ./ scale).^2 / 2); k = k / sum(k);

%% blurring the image
tic
imablur = conv2(conv2(ima3,k,'same'),k','same');
t=toc; 
fprintf(1,'Image blurring time:         %g sec.\n',t )


figure(3);
image(imablur);
colormap(gray(256));
title('Blurred input image');
axis image
drawnow;


%% Computing gradients and orientations
tic
[dx,dy] = gradient(imablur); 
gradmag = sqrt(dx.^2 + dy.^2);
sintheta = dx ./ gradmag;
costheta = - dy ./ gradmag;
t=toc; fprintf(1,'Gradients and angles computation time: %g sec.\n',t )


figure(4);
imagesc(gradmag);
colorbar;
title('Gradient magnitude at every pixel');
axis image
drawnow;


%% Fitting ellypses
tic
[dist,fit_err] = subpixeledges(gradmag,sintheta,costheta);
t=toc;
fprintf(1,'Ellipses fitting time: %g sec.\n',t )



%% calculate appropriate threshold
tic
gradsort = sort(gradmag(:));
maxgrad = gradsort(end);
thresh = gradsort(round(m*n*hithresh));
t=toc; 
fprintf(1,'Gradient threshold is %g (max grad is %g) computation time: %g sec.\n',thresh, maxgrad, t )


% Edge locations (not at subpixel accuracy yet)
edges = sparse((abs(dist) < 0.501) .* (gradmag > thresh) .* (fit_err < fitthresh));


figure(5);
spy(edges);
title('Edge locations (not subpixel)');
zoom on;
axis image
drawnow;


tic
sintheta = sparse( sintheta .* edges );
costheta = sparse( costheta .* edges );
dist = sparse(dist .* edges);
[posx,posy] = meshgrid(1:n,1:m);
posx = sparse(posx .* edges );
posy = sparse(posy .* edges );

posx = posx  + ( dist .* sintheta);
posy = posy - ( dist .* costheta);
t=toc; fprintf(1,'Construction of output arrays time: %g sec.\n',t )


figure(6);
image(ima);
title('Computed subpixel egdes');
hold on;
quiver(posx,posy,gradmag.*sintheta/20,-gradmag.*costheta/20,0);
hold off;
zoom on;
drawnow;
