掌握 canny 边缘检测原理，完成图像的边缘检测实验，展示每个环节的处理结果（梯度图，NMS，边缘链接）

须满足条件：抑制噪声；精确定位边缘。从数学上表达了三个准则[信噪比准则（低错误率）、定位精度准则、单边缘响应准则]，并寻找表达式的最佳 l 解。属于先平滑后求导的方法。算法基本步骤

使用高斯滤波平滑图像

令 f ( x , y ) 表示数据（输入源数据）， G ( x , y ) 表示二维高斯函数（卷积操作数）， f s ( x , y ) 令 f(x,y)表示数据（输入源数据）， G(x,y)表示二维高斯函数（卷积操作数），f_s(x,y)令 f(x,y)表示数据（输入源数据），G(x,y)表示二维高斯函数（卷积操作数），fs(x,y)为卷积平滑后的图像。

Guess 过程

用坐标点( x , y ) 表示一个 3 × 3 的邻域，设中心点的坐标为 ( 0 , 0 ) (x,y)表示一个3\times3的邻域，设中心点的坐标为(0,0)(x,y)表示一个3×3 的邻域，设中心点的坐标为(0,0)，相邻的点以此类推。计算权重矩阵。设定方差σ 2 = 0.64 \sigma^2=0.64σ2=0.64 的值，将对应各个坐标点( x , y ) (x,y)(x,y)带入二维高斯公式 G ( x , y ) G(x,y)G(x,y)中，得到一个权重矩阵，归一化权重矩阵(矩阵中各个点除以权重之和)，得到标准的权重矩阵，即高斯模板。

计算高斯模糊。设在一幅图像中的 3 × 3 3\times33×3 区域内，用各像素点的灰度值乘以对应点的权重。

将得到的 9 个值求和，就是中心点的高斯模糊值。

简单来说就是 Guess 模板在原始图像中进行移位、相乘、相加的过程。

计算幅值图像

角度图像求变化率时，对于一元函数，即求导；对于二元函数，求偏导。数字图像处理中，用一阶有限差分近似求取灰度值的梯度值(变化率)。求灰度的变化率，分别取 x 和 y 方向上相邻像素做差，代替求取 x 和 y 方向一阶偏导)

对幅值图像进行非极大值抑制

步骤： 1. 对角度在四个方向进行划分，形成新的角度图； 2. 根据角度图(表征着梯度的方向)，对幅值进行非极大值抑制。首先将角度划分成四个方向范围:水平( 0 ° ) 、− 45 ° 、垂直 ( 90 ° ) 、+ 45 ° (0°)、-45°、垂直(90°)、+45°(0°)、−45°、垂直(90°)、+45°。接着讨论对 3×3 区域的四个基本边缘方向进行非极大值抑制。做法：若中心点（即：访问点）在沿其方向上邻域的梯度幅值最大，则保留；否则，抑制。

双阈值检测和连接边缘

选取高阈值 T_H 和低阈值 T_L，比率为 2:1 或 3:1。（一般取 T_H=0.3 或 0.2,T_L=0.1）取出非极大值抑制后的图像中的最大梯度幅值，重新定义高低阈值。将小于 T_L 的点抛弃赋 0；将大于 T_H 的点立即标记（这些点就是边缘点），赋 1
将大于 T_L，小于 T_H 的点使用 8 连通区域确定（即：只有与 T_H 像素连接时才会被接受，成为边缘点，赋 1）注意：双阈值做法是将候选像素点拼接成轮廓，轮廓的形成时对这些像素运用滞后性阈值。

clear all;
close all;
clc;

img=imread('testc.jpg');
imshow(img);
[m n]=size(img);
img=double(img);

%%canny边缘检测的前两步相对不复杂，所以我就直接调用系统函数了
%%高斯滤波
w=fspecial('gaussian',[5 5]);
img=imfilter(img,w,'replicate');
figure;
imshow(uint8(img))

%%sobel边缘检测
w=fspecial('sobel');
img_w=imfilter(img,w,'replicate');      %求横边缘
w=w';
img_h=imfilter(img,w,'replicate');      %求竖边缘
img=sqrt(img_w.^2+img_h.^2);        %注意这里不是简单的求平均，而是平方和在开方。我曾经好长一段时间都搞错了
figure;
imshow(uint8(img))

%%下面是非极大抑制
new_edge=zeros(m,n);
for i=2:m-1
    for j=2:n-1
        Mx=img_w(i,j);
        My=img_h(i,j);
        
        if My~=0
            o=atan(Mx/My);      %边缘的法线弧度
        elseif My==0 && Mx>0
            o=pi/2;
        else
            o=-pi/2;            
        end
        
        %Mx处用My和img进行插值
        adds=get_coords(o);      %边缘像素法线一侧求得的两点坐标，插值需要       
        M1=My*img(i+adds(2),j+adds(1))+(Mx-My)*img(i+adds(4),j+adds(3));   %插值后得到的像素，用此像素和当前像素比较 
        adds=get_coords(o+pi);%边缘法线另一侧求得的两点坐标，插值需要
        M2=My*img(i+adds(2),j+adds(1))+(Mx-My)*img(i+adds(4),j+adds(3));   %另一侧插值得到的像素，同样和当前像素比较
        
        isbigger=(Mx*img(i,j)>M1)*(Mx*img(i,j)>=M2)+(Mx*img(i,j)<M1)*(Mx*img(i,j)<=M2); %如果当前点比两边点都大置1
        
        if isbigger
           new_edge(i,j)=img(i,j); 
        end        
    end
end
figure;
imshow(uint8(new_edge))

%%下面是滞后阈值处理
up=120;     %上阈值
low=100;    %下阈值
set(0,'RecursionLimit',10000);  %设置最大递归深度
for i=1:m
    for j=1:n
      if new_edge(i,j)>up &&new_edge(i,j)~=255  %判断上阈值
            new_edge(i,j)=255;
            new_edge=connect(new_edge,i,j,low);
      end
    end
end
figure;
imshow(new_edge==255)
复制代码

function nedge=connect(nedge,y,x,low)       %种子定位后的连通分析
    neighbour=[-1 -1;-1 0;-1 1;0 -1;0 1;1 -1;1 0;1 1];  %八连通搜寻
    [m n]=size(nedge);
    for k=1:8
        yy=y+neighbour(k,1);
        xx=x+neighbour(k,2);
        if yy>=1 &&yy<=m &&xx>=1 && xx<=n  
            if nedge(yy,xx)>=low && nedge(yy,xx)~=255   %判断下阈值
                nedge(yy,xx)=255;
                nedge=connect(nedge,yy,xx,low);
            end
        end        
    end 

end
复制代码