zj_canny.m

canny算子,很有用的一个程序.容易掌握.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%自己编的Canny算子的程序
%I为输入图像,T1为双阈值处理时较大的阈值，T2为双阈值处理时较小的阈值
%2008.6.26
%author:yuanbaofeng
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=zj_canny(I,T1)
T2=0.5*T1;%T2为双阈值处理时较小的阈值，由T1求得
[m,n]=size(I);
%第一步，首先对输入图像进行加噪
%J=imnoise(I,'salt & pepper',0.02);%对图像加上椒盐噪声
%subplot(231);
%imshow(I);
%subplot(232);
%imshow(J);
%第二步，对加噪图像进行高斯平滑
tic;
%A=[1 1 2 2 2 1 1;1 2 2 4 2 2 1;2 2 4 8 4 2 2;2 4 8 16 8 4 2;2 2 4 8 4 2 2;1 2 2 4 2 2 1;1 1 2 2 2 1 1];
%k1=conv2(I,A,'full');
%k1=mat2gray(k1);
%k1=k1*255;
%k1=uint8(k1);
%k1=k1(4:67,4:67);
k1=filter2(fspecial('gaussian',7),I);%对图像进行3*3模版的高斯平滑
%subplot(233);
%imshow(uint8(k1));
%第三步，对平滑图像求梯度幅度和方向
%k1=I;
%k1=I;
Ex=zeros(m-1,n-1);%初始化x方向差分矩阵
Ey=zeros(m-1,n-1);%初始化y方向差分矩阵
%fu_zh=zeros(m-1,n-1);%初始化梯度幅值矩阵
%fa_x=zeros(m-1,n-1);%初始化梯度方向矩阵
%orient=zeros(m-1,n-1);%初始化方向范围矩阵
global flag;
flag=zeros(m-1,n-1);%初始化一个标记矩阵
%求解x,y方向的差分
for i=1:m-1
    for j=1:n-1
        Ex(i,j)=(k1(i,j+1)-k1(i,j)+k1(i+1,j+1)-k1(i+1,j))/2;
        Ey(i,j)=(k1(i,j)-k1(i+1,j)+k1(i,j+1)-k1(i+1,j+1))/2;
    end
end

%求解幅值和方向
Ex1=Ex(:);
Ey1=Ey(:);
% for i=1:m-1
%     for j=1:n-1
%         fu_zh(i,j)=abs(Ex(i,j))+abs(Ey(i,j));
%         fa_x(i,j)=atan(Ey(i,j)/(Ex(i,j))+eps);
%     end
% end
for i=1:(m-1)*(n-1)
      fu_zh(i)=abs(Ex1(i))+abs(Ey1(i));
      fa_x(i)=atan(Ey(i)/(Ex(i))+eps);  
end
fu_zh=reshape(fu_zh,m-1,n-1);
%fa_x=reshape(fa_x,m-1,n-1);
%非极大值抑制法
for i=1:(m-1)*(n-1)
    %for j=1:n-1
        %switch fa_x(i)
            if fa_x(i)>-0.125*pi&fa_x(i)<0.125*pi%求梯度方向范围
               orient1(i)=1;continue;
            end
            if fa_x(i)>-0.375*pi&fa_x(i)<-0.125*pi
               orient1(i)=4;continue;
            end
            if fa_x(i)>-0.5*pi&fa_x(i)<-0.375*pi|fa_x(i)>0.375*pi&fa_x(i)<0.5*pi
               orient1(i)=3;continue;
            end
            if fa_x(i)>0.125*pi&fa_x(i)<0.375*pi
               orient1(i)=2;continue;
            else  
               orient1(i)=3;continue;
        
            end

end
orient=reshape(orient1,m-1,n-1);
for i=2:m-2
    for j=2:n-2
        switch orient(i,j)
            case 1
                if fu_zh(i,j)<fu_zh(i,j-1)|fu_zh(i,j)<fu_zh(i,j+1)%和梯度方向的的像素比较大小，判零或不变
                    fu_zh(i,j)=0;
                else
                    ;
                end
            case 2
                if fu_zh(i,j)<fu_zh(i+1,j-1)|fu_zh(i,j)<fu_zh(i-1,j+1)
                    fu_zh(i,j)=0;
                else
                    ;
                end
            case 3
                if fu_zh(i,j)<fu_zh(i-1,j)|fu_zh(i,j)<fu_zh(i+1,j)
                    fu_zh(i,j)=0;
                else
                    ;
                end
            case 4
                if fu_zh(i,j)<fu_zh(i-1,j-1)|fu_zh(i,j)<fu_zh(i+1,j+1)
                    fu_zh(i,j)=0;
                else
                    ;
                end
        end
    end
end
fu_zh(1,:)=0;
fu_zh(:,1)=0;
fu_zh(m-1,:)=0;
fu_zh(:,n-1)=0;
%非极大值抑制结束
%双阈值法处理
%找两个阈值，一个大T1一个小T2，T2=0.4*T1
fu_zh1=fu_zh(:);   %求得大阈值处理的图像
for i=1:length(fu_zh1)
    if fu_zh1(i)>T1
        fu_zh1(i)=255;
    else 
        fu_zh1(i)=0;
    end
end
fu_zh1=reshape(fu_zh1,m-1,n-1);%重构矩阵
fu_zh2=fu_zh(:);    %求得小阈值处理的图像
for i=1:length(fu_zh2)
    if fu_zh2(i)>T2
        fu_zh2(i)=255;
    else
        fu_zh2(i)=0;
    end
end
fu_zh2=reshape(fu_zh2,m-1,n-1);%重构矩阵
%跟踪边缘连接
for i=1:m-1
    for j=1:n-1
    if fu_zh1(i,j)~=0&flag(i,j)~=1
        flag(i,j)=1;      %如果为非零元素，标记为1
        r=i;
        v=j;
        R=1;
        %R=sum(fu_zh1(i-1:i+1,j-1:j+1));
        while R|(fu_zh1(r,v)~=0&flag(r,v)==0)
        [k,l]=gen_zh1(fu_zh1,i,j);%跟踪轮廓
        [r,v]=gen_zh2(fu_zh2,k,l); %跟踪图像2
        %gen_zh1(fu_zh1,r,v);%返回图像1中跟踪
        R=0;
        end
    end
    end
end
toc;
flag=mat2gray(flag);
imshow(flag);