图像处理黄金算子：Canny算子手写体！！

首先写算法 推荐 听 周杰伦 的 歌，本人 亲身 体验，周杰伦 的 歌会提供最多  的 灵感。

---

其实不是闲，因为老师让写。

还是只讲干货，但这次由于算子原理很牛B，讲下原理和我的理解。

首先，高斯滤波直接 用FS函数生成核，直接在时域滤波就OK了，因为整个程序下来都不涉及频域的计算，所以还是比较顺利，改了半个小时的bug就跑好了。

然后，使用sobel求梯度（幅值和角度，这是两张图）。

非极大值抑制，我使用的是“冈萨雷斯”里的方法，直接用的22.5度系列值，这样可能比插值要准，因为我的直觉里插值在求边缘效果不会很好。 抑制完以后放在一个新矩阵里或者用原来的矩阵，注意所有的抑制后的值都要用，在这里所有的值都要被用来做双阈值处理。非极大值抑制后的矩阵暂且记为A

双阈值：把所有的A里的值分为三段，两个阈值来控制，此时再准备一个矩阵记为mid，mid里位置为大阈值以上的值置1，小阈值一下的置0，中间的在下面处理。

抑制孤立低阈值点：八邻域内有强边缘值，则若边缘置1（注意一定要在mid里置1），千万不要在A里修改像素值！！

mymain：（iCanny）注释真的超级良心的

function imgCannyout=iCanny(I,topvalue,bottomvalue,m,n,sigma)
%%Gauss平滑滤波器
%Gauss卷积核
gausskernel=fspecial('gaussian',[m,n],sigma);         %由Fspecial生成高斯核
%图像
IG=rgb2gray(I);                                      %灰度图直接注释掉
[M,N]=size(IG);
%imgnew=imgplus(IG,(m-1)/2,(n-1)/2);
imgnew=IG;
%Gauss滤波
imgData=zeros(M,N);
for ii=2:M-1
    for jj=2:N-1
        window=imgnew(ii-(m-1)/2:ii+(m-1)/2,jj-(m-1)/2:jj+(m-1)/2);   %取窗口
        Gaussfilter=double(window).*gausskernel;
        imgData(ii,jj)=sum(Gaussfilter(:));                           %卷积
    end
end
%计算梯度值，并生成梯度矩阵（使用sobel算子）
sx=[-1,0,1;-2,0,2;-1,0,1];
sy=[1,2,1;0,0,0;-1,-2,-1];
[m,n]=size(sx);
%imgnew_=imgplus(imgData,(m-1)/2,(n-1)/2);
imgnew_=imgData;
gx=zeros(M,N);
gy=zeros(M,N);
gm=zeros(M,N);
direction=zeros(M,N);
%计算梯度值
for iii=2:M-1
    for jjj=2:N-1
        window_=imgnew_(iii-(m-1)/2:iii+(m-1)/2,jjj-(n-1)/2:jjj+(n-1)/2);
        dx=double(window_).*sx;
        dy=double(window_).*sy;
       % dx=dx';
        dx=sum(dx(:));
        dy=sum(dy(:));
        gx(iii,jjj)=dx;
        gy(iii,jjj)=dy;
        gm(iii,jjj)=sqrt(dx^2+dy^2);                   %梯度幅值
        %计算梯度值
        angle=atand(dy/dx);                           %梯度角度值，按照书上角度分界取值
        if(angle>=67.5 && angle<90)
            direction(iii,jjj)=2;
        elseif(angle>=-67.5 && angle<-22.5)
            direction(iii,jjj)=3;
        elseif(angle>=-22.5 && angle<22.5)
            direction(iii,jjj)=0;
        elseif(angle>=22.5 && angle<67.5)
            direction(iii,jjj)=1;
        else
            direction(iii,jjj)=2;
        end
    end
end
%%
%非极大值抑制
gradmthen=zeros(M,N);
mid=zeros(M,N);
for iiii=2:M-1
    for jjjj=2:N-1
        %八个方向列好避免出错
        t=gm(iiii-1,jjjj);                     
        b=gm(iiii+1,jjjj);
        r=gm(iiii,jjjj+1);
        l=gm(iiii,jjjj-1);
        tr=gm(iiii-1,jjjj+1);
        tl=gm(iiii-1,jjjj-1);
        lb=gm(iiii+1,jjjj-1);
        rb=gm(iiii+1,jjjj+1);
            %把要比较的两个值取出来
         if direction(iiii,jjjj) == 0
            gradm1 = r;
            gradm2 = l;
         elseif direction(iiii,jjjj) == 1
            gradm1 = tr;
            gradm2 = lb;
          elseif direction(iiii,jjjj) == 2
            gradm1 = b;
            gradm2 = t;
         else direction(iiii,jjjj) == 3
            gradm1 = tl;
            gradm2 = rb;
         end
         %抑制
         if gm(iiii,jjjj)>=gradm1 && gm(iiii,jjjj)>=gradm2
             gradmthen(iiii,jjjj)=gm(iiii,jjjj);
         else
             gradmthen(iiii,jjjj)=0;
         end
    end
end
%经过抑制之后要进行阈值处理，此时需要（必须）用到mid矩阵了
for i=2:M-1
    for j=2:N-1
        if gradmthen(i,j)>=topvalue;
            gradmthen(i,j)=topvalue;
            mid(i,j)=1;                       %强边缘为1，弱边缘先保留，但是mid里不置1 ！！
        elseif gradmthen(i,j)>=bottomvalue
             gradmthen(i,j)=bottomvalue;
        else
             gradmthen(i,j)=0;
        end
    end
end
        %专门判断弱边缘
  for i_=2:M-1
    for j_=2:N-1
    if gradmthen(i_,j_)==bottomvalue
        window_=gradmthen(i_-1:i_+1,j_-1:j_+1);
        n=length(find(window_==topvalue));          %找出周围八个点中强边缘的个数 
        if n>=1
            mid(i_,j_)=1;
        end
        end
    end
  end
  imgCannyout=mid;
  imshow(imgCannyout);

如果担心边界溢出影响效果，可以用我main里注释掉到的imgplus（）：

如下：

function imnew__ = imgplus(img,m,n)
m=n;
[M,N] = size(img);                       %控制台输入加噪图像'im2'
s=(m-1)/2;
imnew__=zeros(M+s,N+s);                     %构造扩充边缘矩阵
imnew__(s+1:s+M,s+1:s+N)=img;              
for ii=1:s
imnew__(ii,s+1:s+N)=img(1,:);              %上方行
imnew__(ii+M+s,s+1:s+N)=img(M,:);          %下方行
imnew__(s+1:s+M,ii)=img(:,1);              %左部列
imnew__(s+1:s+M,ii+N+s)=img(:,N);          %右部列
end
for iii=1:s
    for jjj=1:s
        imnew__(iii,jjj)=img(1,1);        %左上角
        imnew__(s+M+iii,jjj)=img(M,1);    %左下角
        imnew__(iii,s+N+jjj)=img(1,N);    %右上角
        imnew__(s+M+iii,s+N+jjj)=img(M,N);%右下角
    end
end

原图：

效果图：由于老师给的lena分辨率太低，有一部分没有连接起来而且图拉大以后很糊。