function [edges,posx,posy,costheta,sintheta,dx,dy,gradmag] = canny(ima,scale,hithresh,fitthresh)

%CANNY          Find canny edges in image at scale s and 
%               thresholds hithresh and lothresh
%
%               Use:    [edges,posx,posy,costheta,sintheta] = canny(ima,scale,hithresh,fitthresh)
%
%               The thresholds are supposed to be between 0 and 1.
%               hithresh is optional. Set equal to 0.7 by default
%               fitthresh is optional. Set equal to 0.1 by default
%
%               After running canny to see the edges try command:
%
%               spy(edges) 
%
%               to get a quiver display at subpixel accuracy try:
%
%               quiver(posx,posy,costheta,sintheta)
%
%               (it runs into an error but manages to do the job
%               try also the quiversparse function)
%


% Conventions: walking along an edge the brighter side is at the right
%               and distances are positive when walking away from the
%               the edge on the dark side, while they are negative
%               when walking away on the bright side.
%               Angles increase with leftwards rotation.
%
%               BE CAREFUL: Matlab inverts the y axis and threfore
%               the angles and other conventions look inverted on the
%               screen


%       (c) Copyright Pietro Perona, 
%       California Inst. of Technology
%                July 8 1993
%       modified May 13 1994
%               

if( (fopen('subpixeledges.m','r')) == -1 )
        error('You need a copy of the file subpixeledges.m');
end

if(nargin==2)
        hithresh = 0.7;
        fitthresh = 0.1;
end;
if(nargin==3)
        fitthresh = 0.4;
end;
lothresh = hithresh/2;

format short
format compact

[m,n] = size(ima);

%% 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(ima,k,'same'),k','same');
%t=toc; fprintf(1,'Image blurring time:         %g sec.\n',t )



%% Computing gradients and orientations
%tic
[dx,dy] = gradient(imablur); 
gradmag = sqrt(dx.^2 + dy.^2);

sintheta = dx ./ gradmag;
costheta = - dy ./ gradmag;

%% Theta is commented out to increase speed
%theta1 = atan(sintheta ./ costheta);
%theta = theta1 + (pi .* sign(sintheta)) .* (costheta<0);
%t=toc; fprintf(1,'Gradients and angles computation time: %g sec.\n',t )


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



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


%tic
%image_borders = sparse(
edges = sparse((abs(dist) < 0.501) .* (gradmag > thresh) .* (fit_err < fitthresh));
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 );

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

clear global sizeimage gradmag hithresh;
format loose

