图像分割——边缘检测——Canny算法（Matlab）

clc;
clear all;
close all;
I=im2double(imread('D:\Gray Files\10-26.tif'));
[M,N]=size(I);
%%
%=============================边缘检测(五)=================================
% Canny Edge Detector
%-------------------------用高斯低通滤波器平滑图像-------------------------
%建筑图像所设参数
% n=25;
% sigma=4;
%头颅扫描图像所设参数
n=13;
sigma=2;
%lena测试图像所设参数
% n=5;
% sigma=1;
type='symmetric';
f_s=GaussianBlur(n,I,sigma,type);
% imshow(f_s)
% imwrite(f_s,'D:\Gray Files\lena-test.jpg','jpg');
%-----------------------计算平滑后的图像梯度和角度-------------------------
n_l=1;
%Sobel算子
s_y=[-1 -2 -1;
    0 0 0;
    1 2 1];
s_x=[-1 0 1;
    -2 0 2;
    -1 0 1];
%定义梯度和角度
gx=zeros(M,N);
gy=zeros(M,N);
f_s_pad=padarray(f_s,[n_l,n_l],'replicate');
for i=1:M
    for j=1:N
        Block=f_s_pad(i:i+2*n_l,j:j+2*n_l);
        gx(i,j)=sum(sum(Block.*s_x));
        gy(i,j)=sum(sum(Block.*s_y));        
    end
end
type='replicate';
gx=GaussianBlur(n,gx,sigma,type);
gy=GaussianBlur(n,gy,sigma,type);
M_s=sqrt(gx.^2+gy.^2);
M_s=M_s/max(M_s(:));
a_s=atan2(gy,gx)*180/pi;
% imshow(M_s)
%-----------------------对梯度图像进行非极大值抑制-------------------------
n_l=1;
%定义非极大值抑制图像
g_N=M_s;
M_s_pad=padarray(M_s,[n_l,n_l],'replicate');
for i=1:M
    for j=1:N
        %取出中心点的梯度值
        K=M_s_pad(i+1,j+1);
 
        theta=a_s(i,j);
        if (theta>=0 && theta<=45) ||...
                (theta<-135 && theta>=-180)
            yBot=[M_s_pad(i+1,j+2) M_s_pad(i+2,j+2)];
            yTop=[M_s_pad(i+1,j) M_s_pad(i,j)];
            k=abs(gy(i,j)/M_s_pad(i+1,j+1));
            K1=(yBot(2)-yBot(1))*k+yBot(1);
            K2=(yTop(2)-yTop(1))*k+yTop(1);         
        end    
        if (theta>45 && theta<=90) ||...
                (theta<-90 && theta>=-135)
            yBot=[M_s_pad(i+2,j+1) M_s_pad(i+2,j+2)];
            yTop=[M_s_pad(i,j+1) M_s_pad(i,j)];
            k=abs(gx(i,j)/M_s_pad(i+1,j+1));
            K1=(yBot(2)-yBot(1))*k+yBot(1);
            K2=(yTop(2)-yTop(1))*k+yTop(1);         
        end            
        if (theta>90 && theta<=135) ||...
                (theta<-45 && theta>=-90)
            yBot=[M_s_pad(i+2,j+1) M_s_pad(i+2,j)];
            yTop=[M_s_pad(i,j+1) M_s_pad(i,j+2)];
            k=abs(gx(i,j)/M_s_pad(i+1,j+1));
            K1=(yBot(2)-yBot