Canny算子理解，及Matlab实现

JohnCanny于1986年提出Canny算子，它与Marr（LoG）边缘检测方法类似，也属于是先平滑后求导数的方法。本文对根据上述的边缘检测过程对Canny检测算法的原理进行介绍。 并结合实验，对canny算法进行理解，指出参数选取中容易产生的问题。

canny边缘检测一共四个部分：
1. 用高斯滤波器平滑图像；（图像去噪）
2. 用一阶偏导有限差分计算梯度幅值和方向；（特征增强）
3. 对梯度幅值进行非极大值抑制 ；（边缘检测）
4. 用双阈值算法检测和连接边缘。（形态学处理）

平滑去噪

canny边缘检测的前两步相对不复杂，所以我就直接调用系统函数了。

IM=imread('2.tif');
IM=imnoise(IM); %imnoise对uint8类型加噪为0-255，e对double类型加噪为0-1，
[m,n]=size(IM);
IM=double(IM);
%高斯滤波
w=fspecial('gaussian',[9 9]);
img=imfilter(IM,w,'replicate');

值得注意的是，在选取滤波窗口w时，根据我的实验（无噪声情况实验，若有噪声，也会影响后续的非极大值抑制），滤波窗口w边长为偶数时，将导致得到的梯度只有单像素最大幅值。而对边缘分析可知，应该产生双像素的梯度幅值才符合实际情况，所以滤波窗口w的选取一定要为基数边长。

实例如下，当选偶数边长时，幅值结果如下

当选基数边长时，幅值结果如下

计算幅值

%%sobel边缘检测
w_h=[1,2,1;0,0,0;-1,-2,-1];              %Sobel算子，梯度方向朝上为正
img_h=imfilter(img,w_h,'replicate');       %梯度是竖着的边缘,即横边缘
w_w=[-1,0,1;-2,0,2;-1,0,1];              %Sobel算子，梯度方向朝右为正
img_w=imfilter(img,w_w,'replicate');       %梯度是横着的边缘,,即竖边缘
img_grad=sqrt(img_w.^2+img_h.^2);             %梯度的绝对值
grad_max=max(max(img_grad));
img_grad=img_grad/grad_max; %归一化

阈值选取

阈值用于后续的边缘连接处理，高阈值代表强边缘，低阈值代表弱边缘，在OpenCV中，阈值需要人工选取。在Matlab中，阈值根据梯度幅值的直方图选取，默认选择幅值大小在前70%的最后一个幅值作为高阈值，低阈值与高阈值的比例为4：10.在我的程序中，由于需要检测很强的边缘，故直方图比例设置得很高。

值得一提的是，强边缘的阈值设置后并不需要包括所有边缘，只需要在弱边缘上有“点状”的分布即可，后续会根据与弱边缘的8连通域情况，得到效果很好的边缘。 另外，对得到的强边缘进行开操作和闭操作可以得到很好的去噪效果，而同样的操作却不能对弱边缘使用。

弱边缘图像：

强边缘图像：

代码如下：

k=1.2;
PercentOfPixelsNotEdges = 1-k*(m+n)/(m*n);%0.995; %强边缘的比例
ThresholdRatio = 0.52;  %强弱边缘的比例
thresh=[ ];
if isempty(thresh)%通过直方图自动计算高低阈值大小
    counts=imhist(img_grad, 256);
    highThresh = find(cumsum(counts) > PercentOfPixelsNotEdges*m*n,1,'first') / 256; %PercentOfPixelsNotEdges=0.8，即不是边界的比例
    lowThresh = ThresholdRatio*highThresh;
    thresh = [lowThresh highThresh];
elseif length(thresh)==1
    highThresh = thresh;
    if thresh>=1
        error(message('images:edge:thresholdMustBeLessThanOne'))
    end
    lowThresh = ThresholdRatio*thresh;
    thresh = [lowThresh highThresh];
elseif length(thresh)==2
    lowThresh = thresh(1);
    highThresh = thresh(2);
    if (lowThresh >= highThresh) || (highThresh >= 1)
        error(message('images:edge:thresholdOutOfRange'))
    end
end

非极大值抑制

在John Canny提出的Canny算子的论文中，非最大值抑制就只是在0、90、45、135四个梯度方向上进行的，每个像素点梯度方向按照相近程度用这四个方向来代替。这种情况下，非最大值抑制所比较的相邻两个像素就是：

1） 0：左边 和 右边
2）45：右上 和 左下
3）90： 上边 和 下边
4）135： 左上 和 右下

这样做的好处是简单， 但是这种简化的方法无法达到最好的效果， 因为，自然图像中的边缘梯度方向不一定是沿着这四个方向的。因此，就有很大的必要进行插值，找出在一个像素点上最能吻合其所在梯度方向的两侧的像素值。

非极大值抑制是canny算法的精髓，对此理解不深的同学，可以参考(http://blog.csdn.net/kezunhai/article/details/11620357) 这里的解释，解释得非常到位。

值得一提的是，为了在噪声中得到真实的双像素边缘，可以在非极大值抑制时不选取最大值为1，只要大小为周围最大值的80%即可。

代码如下：

%%下面是非极大值抑制
weak_edge2=zeros(m+2,n+2);
img_grad2 = padarray(img_grad,[1 1],'replicate'); %为了计算边界上的点，要扩充一下
img_h2 = padarray(img_h,[1 1],'replicate');
img_w2 = padarray(img_w,[1 1],'replicate');
for i=2:m+1         % 图像边缘的像素不能计算
    for j=2:n+1

        if img_grad2(i,j)<=lowThresh
            M1=2;
            M2=2;%这时检测点肯定是不大于1的，故isbigger=0

        else %img_grad(i,j)~=0
            Mx=img_w2(i,j);%梯度是横着的
            My=img_h2(i,j);%梯度是竖着的
            if abs(My)>=abs(Mx) %y方向梯度比x方向梯度大，梯度方向“竖着”的情况
                g1=img_grad2(i-1,j);%g1是检测点正上方一行
                g3=img_grad2(i+1,j);%g3是检测点正下方一行
                weight=(abs(Mx)/abs(My));
                if Mx*My>=0 %x,y方向梯度同号，梯度线分布在第1、3象限，当Mx*My==0时，Mx=0,weight=0
                    g2=img_grad2(i-1,j+1);%g2是检测点上右方
                    g4=img_grad2(i+1,j-1);%g4是检测点下左方
                end
                if Mx*My<0 %x,y方向梯度异号，梯度线分布在第2、4象限
                    g2=img_grad2(i-1,j-1);%g2是检测点上左方
                    g4=img_grad2(i+1,j+1);%g4是检测点下右方
                end
                M1=g1*(1-weight)+g2*weight;%M1是上方的插值
                M2=g3*(1-weight)+g4*weight;%M2是下方的插值 
            end
            if  abs(My)<abs(Mx)%x方向梯度比y方向梯度大，梯度方向“横着”的情况
                g1=img_grad2(i,j+1);%g1是检测点正右方
                g3=img_grad2(i,j-1);%g3是检测点正左方
                weight=(abs(My)/abs(Mx));
                if Mx*My>=0 %x,y方向梯度同号，梯度线分布在第1、3象限
                    g2=img_grad2(i-1,j+1);%g2是检测点上右方
                    g4=img_grad2(i+1,j-1);%g4是检测点下左方
                end
                if Mx*My<0 %x,y方向梯度异号，梯度线分布在第2、4象限
                    g2=img_grad2(i+1,j+1);%g2是检测点下右方
                    g4=img_grad2(i-1,j-1);%g4是检测点上左方
                end
                M1=g1*(1-weight)+g2*weight;%M1是上方的插值
                M2=g3*(1-weight)+g4*weight;%M2是下方的插值 
            end
        end
        %此处选择与经典canny不同，考虑到噪声干扰，为了找到“双线边缘”，假定误差不超过0.9的都同为最大值
        %也可根据加入噪声水平确定该值
        isbigger=(img_grad2(i,j)>=(0.8*M1))&&(img_grad2(i,j)>=(0.8*M2)); %如果当前点比两边点都大，则置为白色 
        if isbigger
           weak_edge2(i,j)=255;  
        end        
   end
end
weak_edge=weak_edge2(2:m+1,2:n+1);
figure(1);
imshow(weak_edge)
[rstrong,cstrong] = find(img_grad>highThresh & weak_edge);
strong_edge=zeros(m,n);
for i=1:length(rstrong)
    r=rstrong(i);
    c=cstrong(i);
    strong_edge(r,c)=weak_edge(r,c);
end
se=strel('square',2);
strong_edge=imopen(strong_edge,se);
se=strel('square',3);
strong_edge=imclose(strong_edge,se);
figure(2);
imshow(strong_edge)

形态学处理

最后，以强边缘为种子，在弱边缘中寻找其8连通域，作为边缘即可。

原图：

canny算法检测到的边界：

代码如下
[rstrong,cstrong] = find(strong_edge);
edge = bwselect(weak_edge, cstrong, rstrong, 8);

在Matlab “edge（）”函数中，最后还做了一个形态学瘦化处理，edge= bwmorph(edge, ‘thin’, 1); 但是我觉得得到的边缘已经足够细了，双像素边缘在之后对边缘区域的操作中也会更加有用。