$Canny$边缘检测算子

基本原理

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

算法实现步骤

$1、$使用高斯滤波平滑图像

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

$$G(x,y)=\frac{1}{2\pi\sigma^2}e^\frac{-(x^2+y^2)}{2\sigma^2}$$
$$f_s(x,y)=f(x,y)\ast G(x,y)$$

$Guess过程$

• 用坐标点$(x,y)表示一个3\times3的邻域，设中心点的坐标为(0,0)$，相邻的点以此类推。
  
  
  
  

• 计算权重矩阵。设定方差$\sigma=0.8$的值，将对应各个坐标点$(x,y)$带入二维高斯公式$G(x,y)$中，得到一个权重矩阵，归一化权重矩阵(矩阵中各个点除以权重之和)，得到标准的权重矩阵，即高斯模板。

  

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

  
  
• 将得到的9个值求和，就是中心点的高斯模糊值。
  
  
  
  
  
  具体过程：http://www.ruanyifeng.com/blog/2012/11/gaussian_blur.html

卷积过程

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

$2、$计算梯度的幅值图像,角度图像

补充：求变化率时，对于一元函数，即求导；对于二元函数，求偏导。

数字图像处理中，用一阶有限差分近似求取灰度值的梯度值(变化率)。
(即：使差商$(\Delta{f}/\Delta{x})近似取代微商(\partial{f}/\partial{x})。求灰度的变化率，分别取x和y方向上相邻像素做差，代替求取x和y$方向一阶偏导)

例：计算一点$x方向和y$方向的梯度幅值和方向

上图中显示一段直的边缘线段放大后一部分，每个方块代表一个像素点，用一个方框强调点处边缘的幅值和方向。令灰色像素值为$0，白色像素值为1$.
如图关于一点为中心的$3\times3$邻域，使用$Prewitt$卷积模板进行计算：
$$p_x= \begin{bmatrix} -1&-1&-1\ 0&0&0\ 1&1&1 \end{bmatrix} , p_y= \begin{bmatrix} -1&0&1\ -1&0&1\ -1&0&1 \end{bmatrix}$$
根据$x方向和y$方向的卷积模板，可知，在$3\times3$邻域中从底部一行像素值减去顶部一行的像素，得到$x$方向的偏导数(梯度)；同样，从右边一列像素值减去左边一列的像素，得到$y$方向的偏导数。
$x$方向的梯度：
$$g_x=\partial{f}/\partial{x}=(0-0)+(0-1)+(0-1)=-2$$
$y$方向的梯度：
$$g_y=\partial{f}/\partial{x}=(1-0)+(1-0)+(0-0)=2$$
由此，可以得到该点梯度的幅值和方向：
$$M(x,y)=\sqrt{g_x^2+g_y^2}=2\sqrt{2}$$
$$\alpha(x,y)=arctan[g_y/g_x]=135°$$
如下图表示了中心点的梯度向量、方位角以及边缘方向。（任一点的边缘与梯度向量正交）

注意：

• $g_x,g_y,M(x,y)和\alpha(x,y)都是与原图像大小相同的图像。$
• $\alpha(x,y)是为了确 定梯度方向，以及确定非极大值对比方向；M(x,y)用来确定是否为边缘点，幅值越大，表示像素区 域灰度变化越明显，更能代表边缘点。$

$Canny$算子卷积模板：

$$p_x= \begin{bmatrix} -1&1\ -1&1 \end{bmatrix} , p_y= \begin{bmatrix} 1&1\ 1&1 \end{bmatrix}$$
$$p[i,j]=(f[i,j+1]-f[i,j]+f[i+1,j+1]-f[i+1,j])/2$$
$$Q[i,j]=(f[i,j]-f[i+1,j]+f[i,j+1]-f[i+1,j+1])/2$$
$$M[i,j]=\sqrt{p[i,j]^2+Q[i,j]^2}$$
$$\theta=arctan(Q[i,j]/p[i,j])$$

$3、$对幅值图像进行应用非极大值抑制

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

$4、$双阈值检测和连接边缘

• 选取高阈值$T_H和低阈值T_L，比率为2:1或3:1。（一般取T_H=0.3/0.2,T_L=0.1$）
• 取出非极大值抑制后的图像中的最大梯度幅值，重新定义高低阈值。即：$T_H\times{Max},T_L\times{Max}$。(当然可以自己给定)
• $将小于T_L的点抛弃，赋0；将大于T_H的点立即标记（这些点就是边缘点），赋1。$
• $将大于T_L，小于T_H的点使用8连通区域确定（即：只有与T_H像素连接时才会被接受，成为边缘点，赋1）$。

注意：双阈值做法是将候选像素点拼接成轮廓，轮廓的形成时对这些像素运用滞后性阈值。

算法实现

Matlab代码

1.  clear all;
2.  clc;
3.  I = imread(‘rice.png‘);%读图
4.  % I = rgb2gray(I);%灰度转换
5.  I = double(I);%转化为双精度
6.  [H,W] = size(I);%获取图像大小
7.  
8.  %%  Step1：使用高斯滤波平滑图像
9.  
10. B = [1 2 1;2 4 2;1 2 1];%高斯滤波系数
11. B = 1/16.*B;%高斯滤波模板 方差=0.8
12. A = conv2(I,B,‘same‘);%使用高斯模板进行卷积.计算二维卷积,结果与原图像大小相同 
13. 
14. %%  Step2：计算梯度的幅值图像,角度图像.
15. 
16. %Prewitt梯度模板
17. dx = [-1 0 1;-1 0 1;-1 0 1];%x方向的梯度模板
18. dy = [1 1 1; 0 0 0;-1 -1 -1];%y方向的梯度模板
19. gx = conv2(A,dx,‘same‘);%获取x方向的梯度图像.使用梯度模板进行二维卷积,结果与原图像大小相同
20. gy = conv2(A,dy,‘same‘);%获取y方向的梯度图像.使用梯度模板进行二维卷积,结果与原图像大小相同
21. M = sqrt((gx.^2) + (gy.^2));%获取幅值图像.大小与原图像相等.(.^)表示数组乘方
22. a = atan2(gy,gx);%获取弧度，范围：-pi~pi
23. a = a*180/pi;%将弧度转换为角度，得到角度图像，与原图像大小相等.
24. 
25. %%  Step3：对幅值图像进行应用非极大值抑制
26. 
27. %首先将角度划分成四个方向范围:水平(0°)、-45°、垂直(90°)、+45°
28. for i = 1:H
29.     for j = 1:W
30.         if((a(i,j)>=-22.5) && (a(i,j)<0)||(a(i,j)>=0) && (a(i,j)<22.5) || (a(i,j)<=-157.5) && (a(i,j)>=-180)||(a(i,j)>=157.5)&&(a(i,j)<=180))
31.             a(i,j) = 0;
32.         elseif((a(i,j)>=22.5) && (a(i,j)<67.5) || (a(i,j)<=-112.5) && (a(i,j)>-157.5))
33.             a(i,j) = -45;
34.         elseif((a(i,j)>=67.5) && (a(i,j)<112.5) || (a(i,j)<=-67.5) && (a(i,j)>-112.5))
35.             a(i,j) = 90;
36.         elseif((a(i,j)>=112.5) && (a(i,j)<157.5) || (a(i,j)<=-22.5) && (a(i,j)>-67.5))
37.             a(i,j) = 45;  
38.         end
39.     end
40. end
41. %讨论对3x3区域的四个基本边缘方向进行非极大值抑制.获取非极大值抑制图像
42. Nms = zeros(H,W);%定义一个非极大值图像
43. for i=2:H-1
44.     for j=2:W-1
45.         if (a(i,j)==0 && M(i,j) == max([M(i,j), M(i,j+1), M(i,j-1)]))
46.             Nms(i,j) = M(i,j);
47.         elseif (a(i,j)==-45 && M(i,j) == max([M(i,j), M(i+1,j-1), M(i-1,j+1)]))
48.             Nms(i,j) = M(i,j);
49.         elseif (a(i,j)==90 && M(i,j) == max([M(i,j), M(i+1,j), M(i-1,j)]))
50.             Nms(i,j) = M(i,j);
51.         elseif (a(i,j)==45 && M(i,j) == max([M(i,j), M(i+1,j+1), M(i-1,j-1)]))
52.             Nms(i,j) = M(i,j);
53.         end;
54.     end;
55. end;
56. 
57. %%  Step4:双阈值检测和连接边缘
58. 
59. DT = zeros(H,W);%定义一个双阈值图像
60. TL = 0.1 * max(max(Nms));%低阈值
61. TH = 0.3 * max(max(Nms));%高阈值
62. for i = 1  : H
63.     for j = 1 : W
64.         if (Nms(i, j) < TL)
65.             DT(i,j) = 0;
66.         elseif (Nms(i, j) > TH)
67.             DT(i,j) = 1 ;
68.         %对TL<Nms(i, j)<TH 使用8连通区域确定
69.         elseif ( Nms(i+1,j)>TH || Nms(i-1,j)>TH || Nms(i,j+1)>TH || Nms(i,j-1)>TH || Nms(i-1, j-1)>TH || Nms(i-1, j+1)>TH || Nms(i+1, j+1)>TH || Nms(i+1, j-1)>TH)
70.             DT(i,j) = 1;
71.         end;
72.     end;
73. end;
74. figure, imshow(DT); %最终的边缘检测为二值图像

效果图