function e = canny(image, sigma)
% obtain gradients

fx = gD(image, sigma, 1, 0);
fy = gD(image, sigma, 0, 1);
%gradient norm:
fw = sqrt(fx .^ 2 + fy .^ 2);
%directions:
ftheta = atan2(fy, fx);

%test
%fw = fw - max(max(fw))/2;

%allocate mask
mask = zeros(size(image));

%suppress non-maximums
for a = 2:size(image,1)-1
    for b = 2:size(image,2)-1
        N = fw(a-1:a+1, b-1:b+1);
        mask(a,b) = checkNeighborhood(N, ftheta(a,b));
    end
end

e = mask .* fw;
