demo_canny.m

img = imread ('lena.jpg');
img = rgb2gray(img); 
img = double (img);

% Value for high and low thresholding
threshold_low = 0.035;
threshold_high = 0.175;
 
%% Gaussian filter definition (https://en.wikipedia.org/wiki/Canny_edge_detector)
G = [2, 4, 5, 4, 2; 4, 9, 12, 9, 4;5, 12, 15, 12, 5;4, 9, 12, 9, 4;2, 4, 5, 4, 2];
G = 1/159.* G;
 
%Filter for horizontal and vertical direction
dx = [1 0 -1];
dy = [1; 0; -1];

%% Convolution of image with Gaussian
Gx = conv2(G, dx, 'same');
Gy = conv2(G, dy, 'same');
 
% Convolution of image with Gx and Gy
Ix = conv2(img, Gx, 'same');
Iy = conv2(img, Gy, 'same');


%% Calculate magnitude and angle
magnitude = sqrt(Ix.*Ix+Iy.*Iy);
angle = atan2(Iy, Ix);
 
%% Edge angle conditioning
angle(angle<0) = pi+angle(angle<0);
angle(angle>7*pi/8) = pi-angle(angle>7*pi/8);

% Edge angle discretization into 0, pi/4, pi/2, 3*pi/4
angle(angle>=0&angle<pi/8) = 0;
angle(angle>=pi/8&angle<3*pi/8) = pi/4;
angle(angle>=3*pi/8&angle<5*pi/8) = pi/2;
angle(angle>=5*pi/8&angle<=7*pi/8) = 3*pi/4;

%% initialize the images
[nr, nc] = size(img);
edge = zeros(nr, nc);
 
%% Non-Maximum Supression
edge = non_maximum_suppression(magnitude, angle, edge); 
 
edge = edge.*magnitude;

%% Hysteresis thresholding
% for weak edge
threshold_low = threshold_low * max(edge(:));
% for strong edge
threshold_high = threshold_high * max(edge(:));  
linked_edge = zeros(nr, nc); 
linked_edge = hysteresis_thresholding2(threshold_low, threshold_high, linked_edge, edge);

non_maximum_suppression.m

function edge = non_maximum_suppression(magnitude, angle, edge)
[nr,nc] = size(edge);
for y = 2: nr-1
    for x = 2: nc-1
        switch angle(y,x)
            case 0
                if magnitude(y,x) >= max(magnitude(y,x-1),magnitude(y,x+1))
                    edge(y,x) = 1;
                end
            case pi/4
                if magnitude(y,x) >= max(magnitude(y-1,x-1),magnitude(y+1,x+1))
                    edge(y,x) = 1;
                end
            case pi/2
                if magnitude(y,x) >= max(magnitude(y-1,x),magnitude(y+1,x))
                    edge(y,x) = 1;
                end
            case 3*pi/4
                if magnitude(y,x) >= max(magnitude(y-1,x+1),magnitude(y+1,x-1))
                    edge(y,x) = 1;
                end
        end
    end
end
end

hysteresis_thresholding2.m

function linked_edge = hysteresis_thresholding2(threshold_low, threshold_high, linked_edge, edge)
[nr,nc] = size(edge);
linked_edge = zeros(nr,nc);
queue = zeros(100001,2);
front = 1;
rear = 1;
for y = 2: nr-1
    for x = 2:nc-1
        if edge(y,x) >= threshold_high
            queue(rear,1) = y;
            queue(rear,2) = x;
            rear = rear+1;
            linked_edge(y,x) = 1;
        end
        while front ~= rear
            % pop
            tmp_i = queue(front,1);
            tmp_j = queue(front,2);
            front = front+1;
            % 8邻域寻找弱点并标记linked_edge（加入queue进行BFS）
            if edge(tmp_i - 1,tmp_j - 1) >= threshold_low
                linked_edge(tmp_i - 1, tmp_j - 1) = 1;
                edge(tmp_i - 1,tmp_j - 1) = 0;
                % push
                queue(rear,1) = tmp_i - 1;
                queue(rear,2) = tmp_j - 1;
                rear = rear+1;
            end
            if edge(tmp_i - 1,tmp_j) >= threshold_low
                linked_edge(tmp_i - 1, tmp_j) = 1;
                edge(tmp_i - 1,tmp_j) = 0;
                queue(rear,1) = tmp_i - 1;
                queue(rear,2) = tmp_j;
                rear = rear+1;
            end
            if edge(tmp_i,tmp_j - 1) >= threshold_low
                linked_edge(tmp_i, tmp_j - 1) = 1;
                edge(tmp_i,tmp_j - 1) = 0;
                queue(rear,1) = tmp_i;
                queue(rear,2) = tmp_j - 1;
                rear = rear+1;
            end
            if edge(tmp_i - 1,tmp_j + 1) >= threshold_low
                linked_edge(tmp_i - 1, tmp_j + 1) = 1;
                edge(tmp_i - 1,tmp_j + 1) = 0;
                queue(rear,1) = tmp_i - 1;
                queue(rear,2) = tmp_j + 1;
                rear = rear+1;
            end
            if edge(tmp_i,tmp_j + 1) >= threshold_low
                linked_edge(tmp_i, tmp_j + 1) = 1;
                edge(tmp_i,tmp_j + 1) = 0;
                queue(rear,1) = tmp_i;
                queue(rear,2) = tmp_j + 1;
                rear = rear+1;
            end
            if edge(tmp_i + 1,tmp_j) >= threshold_low
                linked_edge(tmp_i + 1, tmp_j) = 1;
                edge(tmp_i + 1,tmp_j) = 0;
                queue(rear,1) = tmp_i + 1;
                queue(rear,2) = tmp_j;
                rear = rear+1;
            end
            if edge(tmp_i + 1,tmp_j - 1) >= threshold_low
                linked_edge(tmp_i + 1, tmp_j - 1) = 1;
                edge(tmp_i + 1,tmp_j - 1) = 0;
                queue(rear,1) = tmp_i + 1;
                queue(rear,2) = tmp_j - 1;
                rear = rear+1;
            end
            if edge(tmp_i + 1,tmp_j + 1) >= threshold_low
                linked_edge(tmp_i + 1, tmp_j + 1) = 1;
                edge(tmp_i + 1,tmp_j + 1) = 0;
                queue(rear,1) = tmp_i + 1;
                queue(rear,2) = tmp_j + 1;
                rear = rear+1;
            end
        end
    end
end
end