Canny edge dection

edge_detection.m

function output = my_edge(input_image)
    mean_value = mean(mean(input_image));
    [height, width] = size(input_image);
    fft_mid_img = fftshift(fft2(input_image));
    fft_ratio = 12;
    a = round(height/fft_ratio):round((fft_ratio-1)*height/fft_ratio);
    b = round(width/fft_ratio):round((fft_ratio-1)*width/fft_ratio);
%     disp(a);
%     disp(b);
    fft_mid_img(round(height / 2), round(width / 2)) = 0;
    input_image = real(ifft2(ifftshift(fft_mid_img)));
    blurred = gauss_blur(gauss_blur(input_image));
%     blurred = imbilatfilt(input_image);
%     blurred = gauss_blur(medfilt2(input_image));
    subplot(2, 2, 2); imshow(blurred);
    [g, theta] = sobel(blurred);
    
    edge = zeros(height, width);
    for i = 2 : height - 1
        for j = 2 : width - 1
            if (...
                (abs(theta(i, j)) == 90 &&...
                g(i, j) >= g(i+1, j) &&...
                g(i, j) >= g(i-1, j)) ||...
                (theta(i, j) == 0 &&...
                g(i, j) >= g(i, j + 1) &&...
                g(i, j) >= g(i, j - 1)) ||...
                (theta(i, j) == 45 &&...
                g(i, j) >= g(i + 1, j + 1) &&...
                g(i, j) >= g(i - 1, j - 1)) ||...
                (theta(i, j) == -45 &&...
                g(i, j) >= g(i + 1, j - 1) &&...
                g(i, j) >= g(i - 1, j + 1))...
            )
                edge(i, j) = g(i, j);
            else
                edge(i, j) = 0;
            end
        end
    end 
%     upper_thresh = otsu_threshold(input_image);
    upper_thresh = 1.33 * mean_value;
    lower_thresh = 0.66 * mean_value;
%     lower_thresh = upper_thresh / 2;
    disp(upper_thresh);
    disp(lower_thresh);
%     edge(edge > lower_thresh) = 1.0;
    edge((lower_thresh < edge) & (edge < upper_thresh)) = 0.5;
    edge(edge < lower_thresh) = 0;
    edge(edge > upper_thresh) = 1.0;
    edge = filter_low_thresh(edge);
    subplot(2, 2, 3); imshow(edge);
    output = edge;
    
    function filter_ret = filter_low_thresh(filter_img)
        [filter_height, filter_width] = size(filter_img);
        filter_ret = zeros(filter_height, filter_width);
%         filter_ret = filter_img;
        fb = max(filter_height, filter_width) / 20;
        fr = round(max(filter_height, filter_width) / 16);
        for fi = 2:filter_height-1
            for fj = 2:filter_width-1
                if (filter_img(fi, fj) < 0.6 && filter_img(fi, fj) > 0.4)
                    
                    filter_hr = max(1, fi-fr): min(filter_height, fi+fr);
                    filter_wr = max(1, fj-fr): min(filter_width, fj+fr);
%                     disp(filter_wr);
                    filter_area = filter_img(filter_hr, filter_wr);
                    filter_x = sum(sum(filter_area>0.4));
%                     filter_x(ceil(fb/2):filter_height-ceil(fb/2), ...
%                         ceil(fb/2):filter_width-ceil(fb/2)) = 0;
                    if (filter_x > 1)
%                     if (filter_x > 0.6)
                        filter_ret(fi, fj) =1;
                    end 
                elseif (filter_img(fi, fj) > 0.6)
                    filter_ret(fi, fj) = 1;
                end
            end
        end
        
    
    function threshold = otsu_threshold(otsu_img)
        otsu_img = im2uint8(otsu_img) + 1;
        hist = zeros(256, 1);
        [otsu_height, otsu_width] = size(otsu_img);
        for oi = 1:otsu_height
            for oj = 1:otsu_width
                hist(otsu_img(oi, oj)) = hist(otsu_img(oi, oj)) + 1;
            end
        end
        
        weights = 0;
        for oi = 1:256
            weights = weights + hist(oi);
        end
        
        sum = 0.0;
        tmpMax = 0.0;
        threshold = 0;
        background = 0;
        foreground = 0;
        for oi = 1:256
            background = background + hist(oi);
            if (background == 0)
                continue
            end
            foreground = numel(otsu_img) - background;
            if (foreground == 0)
                continue
            end
            sum = sum + oi * hist(oi);
            mean_back = sum / background;
            mean_fore = (weights - sum) / foreground;
            between = background * foreground * (mean_back - mean_fore) ^ 2;
            if (between > tmpMax)
                tmpMax = between;
                threshold = oi;
            end
        end
        threshold = double(threshold - 1) / 255;
        
        
    function [sobel_g, sobel_theta] = sobel(sobel_img)
        sobel_x = [-1 0 1; -2 0 2; -1 0 1];
        sobel_y = [1 2 1; 0 0 0; -1 -2 -1];
        x_ret = imfilter(sobel_img, sobel_x);
        y_ret = imfilter(sobel_img, sobel_y);
        sobel_g = sqrt(x_ret.^2 + y_ret.^2);
        sobel_theta = atan(y_ret./x_ret) * 180 / pi;
        sobel_theta = round(sobel_theta / 45) * 45;
        sobel_theta(isnan(sobel_theta)) = 2;
        
    function gauss_ret = gauss_blur(gauss_img)
        filter_matrix = 1/159 * [
            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];
        gauss_ret = imfilter(gauss_img, filter_matrix);
        
%in this function, you should finish the edge detection utility.
%the input parameter is a matrix of a gray image
%the output parameter is a matrix contains the edge index using 0 and 1
%the entries with 1 in the matrix shows that point is on the edge of the
%image
%you can use different methods to complete the edge detection function
%the better the final result and the more methods you have used, you will get higher scores