计算机视觉（7）————特征检测（边缘／角／斑点）小专题(5)Canny检测

function [imgGray,imgCanny,gradm,jiaodu] = edge_canny(I,gaussDim,sigma,percentOfPixelsNotEdges,thresholdRatio)
%% Gaussian smoothing filter.
m = gaussDim(1);
n = gaussDim(2);
if mod(m,2) == 0 || mod(n,2) == 0
    error('The dimensionality of Gaussian must be odd!');
end
% Generate gaussian convolution kernel.
gaussKernel = fspecial('gaussian', [m,n], sigma);
% Image edge copy.
[m,n] = size(gaussKernel);
[row,col,dim] = size(I);
if dim > 1
    imgGray = rgb2gray(I);
else
    imgGray = I;
end
imgCopy = imgReplicate(imgGray,(m-1)/2,(n-1)/2);
% Gaussian smoothing filter.
imgData = zeros(row,col);
for ii = 1:row
    for jj = 1:col
        window = imgCopy(ii:ii+m-1,jj:jj+n-1);
        GSF = window.*gaussKernel;
        imgData(ii,jj) = sum(GSF(:));
    end
end
%% Calculate the gradient values for each pixel.
% Sobel operator.
dgau2Dx = [-1 0 1;-2 0 2;-1 0 1];
dgau2Dy = [1 2 1;0 0 0;-1 -2 -1];
[m,n] = size(dgau2Dx);
% Image edge copy.
imgCopy = imgReplicate(imgData,(m-1)/2,(n-1)/2);
% To store the gradient and direction information.
gradx = zeros(row,col);
grady = zeros(row,col);
gradm = zeros(row,col);
jiaodu = zeros(row,col);
dir = zeros(row,col); % Direction of gradient.
% Calculate the gradient values for each pixel.
for ii = 1:row
    for jj = 1:col
        window = imgCopy(ii:ii+m-1,jj:jj+n-1);
        dx = window.*dgau2Dx;
        dy = window.*dgau2Dy;
        dx = dx'; % Make the sum more accurate.
        dx = sum(dx(:));
        dy = sum(dy(:));
        gradx(ii,jj) = dx;
        grady(ii,jj) = dy;
        gradm(ii,jj) = sqrt(dx^2 + dy^2);
        % Calculate the angle of the gradient.
        theta = atand(dy/dx) + 90; % 0~180.
        jiaodu(ii,jj)=theta;
        % Determine the direction of the gradient.
        if (theta >= 0 && theta < 45)
            dir(ii,jj) = 2;
        elseif (theta >= 45 && theta < 90)
            dir(ii,jj) = 3;
        elseif (theta >= 90 && theta < 135)
            dir(ii,jj) = 0;
        else
            dir(ii,jj) = 1;
        end
    end
end
% Normalize for threshold selection.
magMax = max(gradm(:));
if magMax ~= 0
    gradm = gradm / magMax;
end
%% Plot 3D gradient graph.
% [xx, yy] = meshgrid(1:col, 1:row);
% figure;
% surf(xx,yy,gradm);
%% Threshold selection.
counts = imhist(gradm, 64);
highThresh = find(cumsum(counts) > percentOfPixelsNotEdges*row*col,1,'first') / 64;
lowThresh = thresholdRatio*highThresh;
%% Non-Maxima Suppression(NMS) Using Linear Interpolation.
gradmCopy = zeros(row,col);
imgBW = zeros(row,col);
for ii = 2:row-1
    for jj = 2:col-1
        E =  gradm(ii,jj+1);
        S =  gradm(ii+1,jj);
        W =  gradm(ii,jj-1);
        N =  gradm(ii-1,jj);
        NE = gradm(ii-1,jj+1);
        NW = gradm(ii-1,jj-1);
        SW = gradm(ii+1,jj-1);
        SE = gradm(ii+1,jj+1);
        % Linear interpolation.
        % dy/dx = tan(theta).
        % dx/dy = tan(90-theta).
        gradValue = gradm(ii,jj);
        if dir(ii,jj) == 0
            d = abs(grady(ii,jj)/gradx(ii,jj));
            gradm1 = E*(1-d) + NE*d;
            gradm2 = W*(1-d) + SW*d;
        elseif dir(ii,jj) == 1
            d = abs(gradx(ii,jj)/grady(ii,jj));
            gradm1 = N*(1-d) + NE*d;
            gradm2 = S*(1-d) + SW*d;
        elseif dir(ii,jj) == 2
            d = abs(gradx(ii,jj)/grady(ii,jj));
            gradm1 = N*(1-d) + NW*d;
            gradm2 = S*(1-d) + SE*d;
        elseif dir(ii,jj) == 3
            d = abs(grady(ii,jj)/gradx(ii,jj));
            gradm1 = W*(1-d) + NW*d;
            gradm2 = E*(1-d) + SE*d;
        else
            gradm1 = highThresh;
            gradm2 = highThresh;
        end
        % Non-Maxima Suppression.
        if gradValue >= gradm1 && gradValue >= gradm2
            if gradValue >= highThresh
                imgBW(ii,jj) = 1;
                gradmCopy(ii,jj) = highThresh;
            elseif gradValue >= lowThresh
                gradmCopy(ii,jj) = lowThresh;
            else
                gradmCopy(ii,jj) = 0;
            end
        else
            gradmCopy(ii,jj) = 0;
        end
    end
end
%% High-Low threshold detection.Double-Threshold.
% If the 8 pixels around the low threshold point have high threshold, then
% the low threshold pixel should be retained.
imgCopy = imgReplicate(imgGray,(m-1)/2,(n-1)/2);

for ii = 3:row-2
    for jj = 3:col-2
        if gradmCopy(ii,jj) == lowThresh
            for x=-2:2
                for y=-2:2
                    if gradmCopy(ii+x,jj+y)==highThresh
                        imgBW(ii,jj) = 1;
                    end
                end
            end                      
        end
    end
end
if magMax ~= 0
    gradm = gradm * magMax;
end
count=1;
c=0;
while ~isempty(count)
    c=c+1;
    num=0;
    disp(length(count))
    count=[];
    disp(length(count))
    for ii=2:row-1
        for jj=2:col-1
            if imgBW(ii,jj)==1 
                for x=-1:1
                    for y=-1:1
                        if (abs(gradm(ii+x,jj+y)-gradm(ii,jj))<5)&&(abs(jiaodu(ii+x,jj+y)-jiaodu(ii,jj))<2)&&imgBW(ii+x,jj+y)==0
                            imgBW(ii+x,jj+y) = 1;
                            num=num+1;
                            count=[count,num,ii+x,ii+y];
                        end
                    end
                end
            end
        end
    end
    disp(count)
end
fprintf( 'c=%d\n ',c)
imgCanny = logical(imgBW);
end
%% Local functions. Image Replicate.
function imgRep = imgReplicate(I,rExt,cExt)
[row,col] = size(I);
imgCopy = zeros(row+2*rExt,col+2*cExt);
% 4 edges and 4 corners pixels.
top = I(1,:);
bottom = I(row,:);
left = I(:,1);
right = I(:,col);
topLeftCorner = I(1,1);
topRightCorner = I(1,col);
bottomLeftCorner = I(row,1);
bottomRightCorner = I(row,col);
% The coordinates of the oroginal image after the expansion in the new graph.
topLeftR = rExt+1;
topLeftC = cExt+1;
bottomLeftR = topLeftR+row-1;
bottomLeftC = topLeftC;
topRightR = topLeftR;
topRightC = topLeftC+col-1;
bottomRightR = topLeftR+row-1;
bottomRightC = topLeftC+col-1;
% Copy original image and 4 edges.
imgCopy(topLeftR:bottomLeftR,topLeftC:topRightC) = I;
imgCopy(1:rExt,topLeftC:topRightC) = repmat(top,[rExt,1]);
imgCopy(bottomLeftR+1:end,bottomLeftC:bottomRightC) = repmat(bottom,[rExt,1]);
imgCopy(topLeftR:bottomLeftR,1:cExt) = repmat(left,[1,cExt]);
imgCopy(topRightR:bottomRightR,topRightC+1:end) = repmat(right,[1,cExt]);
% Copy 4 corners.
for ii = 1:rExt
    for jj = 1:cExt
        imgCopy(ii,jj) = topLeftCorner;
        imgCopy(ii,jj+topRightC) = topRightCorner;
        imgCopy(ii+bottomLeftR,jj) = bottomLeftCorner;
        imgCopy(ii+bottomRightR,jj+bottomRightC) = bottomRightCorner;
    end
end
imgRep = imgCopy;
end
%% End of file.

close all;
clear;
clc;
%%
[filename, pathname] = uigetfile({'*.jpg'; '*.bmp'; '*.gif';'*.png'}, '打开图片');
img= imread([pathname, filename]);
% img = imnoise(img,'salt & pepper',0.01);
%%
[~,~,dim] = size(img);
if dim > 1
    imgCanny1 = edge(rgb2gray(img),'canny'); % system function.
else
    imgCanny1 = edge(img,'canny'); % system function.
end
[imgGray,imgCanny2,gradm,jiaodu] = edge_canny(img,[5,5],1.4,0.7,0.4); % my function.
%%
figure;
imshow(img);
title('original');
%%
figure;
imshow(imgCanny1);
title('system function');
%%
figure;
imshow(imgCanny2);
title('my function');