function [sFinal,thresh] = canny(img, mLow, mHigh, sig)
% Canny edge detector 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function  [sFinal, thresh] =canny(img, mLow, mHigh, sig)
% Applies the canny edge detection algo to the given image
%  img  : the given image (matrix) in color or B/W   
%  mLow : the threshold is set adaptively: low_threshold = mLow* mean_intensity(im_gradient)
%  mHigh: the threshol is set adaptively: high_threshold= mHigh*low_threshold
%  sig  : the value of sigma for the derivative of gaussian operator
%
% The default values for (sig, mLow, mHigh) are (1, 0.5, 2.5)
% The function displays the image and also returns:
%  sFinal       : the final (B/W) image with edges
%  thresh       :  =[lowT, highT] the actual low and high thresholds used
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CS223b HW# 1a: Canny Edge Detector
% Rohit Singh & Mitul Saha:  {rohitsi, mitul}@stanford.edu
%
% Good parameters: 
%     elephant.jpg :  canny(eim,  1.5, 2.6, 1)
%     macbeth.jpg  :  canny(mim, .1, 9, 2); recognises 21 of 24 squares
%     hecface.jpg  :  canny(him, .6, 3, 1); or canny(him, .4,3.5,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (nargin < 1)
  error(' Need a NxNx3 or NxN image matrix');
elseif (nargin ==1)
  mLow = 0.5; mHigh = 2.5; sig = 1;
elseif (nargin == 2)
  mHigh = 2.5; sig = 1;
elseif (nargin == 3)
  sig = 1;
end


%sig = 1;
%mLow = 0.5; 
%mHigh = 2.5;

%img  = imread('macbeth.jpg');
%img  = imread('elephant.jpg');
%img  =  imread('hecface.jpg');

origImage = img;

if (ndims(img)==3)
  img =double(rgb2gray(img));
end

%smoothen ??
%G = gauss(sig);
%img = conv2(img, G'*G,'same');

%CONVOLUTION WITH DERIVATIVE OF GAUSSIAN

dG=dgauss(sig);

[dummy, filterLen] = size(dG);
offset = (filterLen-1)/2;

sy = conv2(img, dG ,'same');
sx = conv2(img, dG','same');

[m, n]=size(img);

% crop off the boundary parts...the places where the convolution was partial
sx = sx(offset+1:m-offset, offset+1:n-offset); 
sy = sy(offset+1:m-offset, offset+1:n-offset); 

% norm of gradient
sNorm = sqrt( sx.^2 + sy.^2 );

% direction of gradient
sAngle = atan2( sy, sx) * (180.0/pi);

% handle divide by zero....
sx(sx==0) = 1e-10;
sSlope = abs(sy ./ sx);

sAorig = sAngle;

%for us, x and x-pi are the same....
y = sAngle < 0;
sAngle = sAngle + 180*y;

% bin the angles into 4 principal directions
% 0-45 45-90 90-135 135-180

binDist =    [-inf 45 90 135 inf];

[dummy, b] = histc(sAngle,binDist);

sDiscreteAngles = b;
[m,n] = size(sDiscreteAngles);

% each pixel is set to either 1,2,3 or 4
% set the boundary pixels to 0, so we don't count them in analysis...
sDiscreteAngles(1,:) = 0;
sDiscreteAngles(end,:)=0;
sDiscreteAngles(:,1) = 0;
sDiscreteAngles(:,end) = 0;

sEdgepoints = zeros(m,n);

sFinal = sEdgepoints;

lowT  = mLow * mean(sNorm(:));
highT = mHigh * lowT;

thresh = [ lowT highT];

% NON-MAXIMAL SUPPRESSION

% for each kind of direction, interpolate.... and see if the current point is
% is the local maximum in that direction

%in the following comments, assume that we are currently at (0,0) and
%we have to interpolate from the 8 surrounding pixels, also, we are
%using MATLAB's single index feature...

%gradient direction: 0-45 i.e. gradDir =1
gradDir = 1;
indxs = find(sDiscreteAngles == gradDir);
slp = sSlope(indxs);

% interpolate between (1,1) and (1,0)
% gDiff1 = Gy/Gx*(magtd(0,0) - magtd(1,1)) + (1- Gy/Gx)*(magtd(0,0)-magtd(1,0))
gDiff1 = slp.*(sNorm(indxs)-sNorm(indxs+m+1)) + (1-slp).*(sNorm(indxs)-sNorm(indxs+1));

% interpolate between (-1,-1) and (-1,0)
% gDiff2 = Gy/Gx*(magtd(0,0) - magtd(-1,-1)) +  (1- Gy/Gx)*(magtd(0,0)-magtd(-1,0))
gDiff2 = slp.*(sNorm(indxs)-sNorm(indxs-m-1)) + (1-slp).*(sNorm(indxs)-sNorm(indxs-1));

okIndxs = indxs( gDiff1 >=0 & gDiff2 >= 0);
sEdgepoints(okIndxs) = 1;


%gradient direction: 45-90 i.e. gradDir =2
gradDir = 2;
indxs = find(sDiscreteAngles == gradDir);
invSlp = 1 ./ sSlope(indxs);

% interpolate between (1,1) and (0,1)
% gDiff1 = (Gx/Gy)*(magtd(0,0) - magtd(1,1)) + (1- Gx/Gy)*(magtd(0,0)-magtd(0,1))
gDiff1 =   invSlp.*(sNorm(indxs)-sNorm(indxs+m+1)) + (1-invSlp).*(sNorm(indxs)-sNorm(indxs+m));

% interpolate between (-1,-1) and (0,-1)
% gDiff2 = (Gx/Gy)*(magtd(0,0) - magtd(-1,-1)) + (1- Gx/Gy)*(magtd(0,0)-magtd(0,-1))
gDiff2 =   invSlp.*(sNorm(indxs)-sNorm(indxs-m-1)) + (1-invSlp).*(sNorm(indxs)-sNorm(indxs-m));

okIndxs = indxs( gDiff1 >=0 & gDiff2 >= 0);
sEdgepoints(okIndxs) = 1;



%gradient direction: 90-135 i.e. gradDir =3
gradDir = 3;
indxs = find(sDiscreteAngles == gradDir);
invSlp = 1 ./ sSlope(indxs);

% interpolate between (-1,1) and (0,1)
% gDiff1 = (Gx/Gy)*(magtd(0,0) - magtd(-1,1)) + (1- Gx/Gy)*(magtd(0,0)-magtd(0,1))
gDiff1 =   invSlp.*(sNorm(indxs)-sNorm(indxs+m-1)) + (1-invSlp).*(sNorm(indxs)-sNorm(indxs+m));

% interpolate between (1,-1) and (0,-1)
% gDiff2 = (Gx/Gy)*(magtd(0,0) - magtd(1,-1)) + (1- Gx/Gy)*(magtd(0,0)-magtd(0,-1))
gDiff2 =   invSlp.*(sNorm(indxs)-sNorm(indxs-m+1)) + (1-invSlp).*(sNorm(indxs)-sNorm(indxs-m));

okIndxs = indxs( gDiff1 >=0 & gDiff2 >= 0);
sEdgepoints(okIndxs) = 1;



%gradient direction: 135-180 i.e. gradDir =4
gradDir = 4;
indxs = find(sDiscreteAngles == gradDir);
slp = sSlope(indxs);

% interpolate between (-1,1) and (-1,0)
% gDiff1 = Gy/Gx*(magtd(0,0) - magtd(-1,1)) + (1- Gy/Gx)*(magtd(0,0)-magtd(-1,0))
gDiff1 = slp.*(sNorm(indxs)-sNorm(indxs+m-1)) + (1-slp).*(sNorm(indxs)-sNorm(indxs-1));

% interpolate between (1,-1) and (1,0)
% gDiff2 = Gy/Gx*(magtd(0,0)-magtd(1,-1)) +  (1- Gy/Gx)*(magtd(0,0)-magtd(1,0))
gDiff2 = slp.*(sNorm(indxs)-sNorm(indxs-m+1)) + (1-slp).*(sNorm(indxs)-sNorm(indxs+1));

okIndxs = indxs( gDiff1 >=0 & gDiff2 >= 0);
sEdgepoints(okIndxs) = 1;


%HYSTERESIS PART...

sEdgepoints = sEdgepoints*0.6;
x = find(sEdgepoints > 0 & sNorm < lowT);
sEdgepoints(x)=0;
x = find(sEdgepoints > 0 & sNorm  >= highT);
sEdgepoints(x)=1;

%sFinal(sEdgepoints>0)=1;

%at this point, if 
%    sNorm(pixel) > lowT then sEdgepoints(pixel)=0
%    highT > sNorm(pixel) > lowT then sEdgepoints(pixel)=0.6
%    sNorm(pixel) > highT then sEdgepoints(pixel)=1

% the idea is this: along the neighbouring pixels in that direction
% add 0.4...so all points that were 0.6 will become 1.0
% see if the number of 1's has changed...
% keep doing while number of 1's doesn't change

oldx = [];
x = find(sEdgepoints==1);
while (size(oldx,1) ~= size(x,1))
  oldx = x;
  v = [x+m+1, x+m, x+m-1, x-1, x-m-1, x-m, x-m+1, x+1];
  sEdgepoints(v) = 0.4 + sEdgepoints(v);
  y = find(sEdgepoints==0.4);
  sEdgepoints(y) = 0;
  y = find(sEdgepoints>=1);
  sEdgepoints(y)=1;
  x = find(sEdgepoints==1);
end
		   
x = find(sEdgepoints==1);

sFinal(x)=1;

figure(1);
imagesc(sFinal); colormap(gray); axis image;