Canny边缘检测
声明：阅读本文需要了解线性代数里面的点乘（图像卷积的原理），高等数学里的二元函数的梯度，极大值定义，了解概率论里的二维高斯分布
1.canny边缘检测原理和简介
2.实现步骤
3.总结
一、 Canny边缘检测算法的发展历史
Canny算子是28岁的John Canny在1986年提出的，该文章发表在PAMI顶级期刊(1986. A computational approach to edge detection. IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 8, 1986, pp. 679-698)。现在老大爷目前（61岁）在加州伯克利做machine learning,主页(http://www.cs.berkeley.edu/~jfc/)，大爷就是大爷。
边缘检测是从图像中提取有用的结构信息的一种技术，如果学过信息论就会知道，一面充满花纹的墙要比一面白墙的信息量大很多，没学过也没关系，直观上也能理解：充满花纹的图像要比单色图像信息更丰富。为什么要检测边缘?因为我们需要计算机自动的提取图像的底层（纹理等）或者高层（时间地点人物等）的信息，边缘可以说是最直观、最容易发现的一种信息了。Canny提出了一个对于边缘检测算法的评价标准，包括：
1) 以低的错误率检测边缘，也即意味着需要尽可能准确的捕获图像中尽可能多的边缘。
2) 检测到的边缘应精确定位在真实边缘的中心。
3) 图像中给定的边缘应只被标记一次，并且在可能的情况下，图像的噪声不应产生假的边缘。
简单来说就是，检测算法要做到：边缘要全，位置要准，抵抗噪声的能力要强。
接下来介绍最经典的canny边缘检测算法，很多边缘检测算法都是在此基础上进行改进的，学习它有利于一通百通。
二、实现步骤
step1：高斯平滑滤波
没有哪张图片是没有噪声的。————鲁迅
滤波是为了去除噪声，选用高斯滤波也是因为在众多噪声滤波器中，高斯表现最好（表现怎么定义的？最好好到什么程度？），你也可以试试其他滤波器如均值滤波、中值滤波等等。一个大小为(2k+1)x(2k+1)的高斯滤波器核（核一般都是奇数尺寸的）的生成方程式由下式给出：
‘
下面是一个sigma = 1.4，尺寸为3x3的高斯卷积核的例子，注意矩阵求和值为1（归一化）：
举个例子：若图像中一个3x3的窗口为A，要滤波的像素点为e，则经过高斯滤波之后，像素点e的亮度值为：
其中*为卷积符号，sum表示矩阵中所有元素相加求和，简单说，就是滤波后的每个像素值=其原像素中心值及其相邻像素的加权求和。图像卷积是图像处理中非常重要且广泛使用的操作，一定要理解好。
其中高斯卷积核的大小将影响Canny检测器的性能。尺寸越大，去噪能力越强，因此噪声越少，但图片越模糊，canny检测算法抗噪声能力越强，但模糊的副作用也会导致定位精度不高，一般情况下，推荐尺寸5*5，3*3也行。
step2： 计算梯度强度和方向
边缘的最重要的特征是灰度值剧烈变化，如果把灰度值看成二元函数值，那么灰度值的变化可以用二元函数的”导数“（或者称为梯度）来描述。由于图像是离散数据，导数可以用差分值来表示，差分在实际工程中就是灰度差，说人话就是两个像素的差值。一个像素点有8邻域，那么分上下左右斜对角，因此Canny算法使用四个算子来检测图像中的水平、垂直和对角边缘。算子是以图像卷积的形式来计算梯度，比如Roberts，Prewitt，Sobel等，这里选用Sobel算子来计算二维图像在x轴和y轴的差分值（这些数字的由来？），将下面两个模板与原图进行卷积，得出x和y轴的差分值图，最后计算该点的梯度G和方向θ
计算梯度的模和方向属于高等数学部分的内容，如果不理解应该补习一下数学基本功，图像处理经常会用到这个概念。
这部分我实现了下，首先了解opencv的二维滤波函数：dst=cv.filter2D(src, ddepth, kernel[, dst[, anchor[, delta[, borderType]]]])
dst： 输出图片
src： 输入图片
ddepth： 输出图片的深度, 详见 combinations,如果填-1，那么就跟与输入图片的格式相同。
kernel： 单通道、浮点精度的卷积核。
以下是默认参数：
anchor：内核的基准点(anchor)，其默认值为(-1,-1)表示位于kernel的中心位置。基准点即kernel中与进行处理的像素点重合的点。举个例子就是在上面的step1中，e=H*A得到的e是放在原像素的3*3的哪一个位置，一般来说都是放在中间位置，设置成默认值就好。
delta ：在储存目标图像前可选的添加到像素的值，默认值为0。（没用过）
borderType：像素向外逼近的方法，默认值是BORDER_DEFAULT,即对全部边界进行计算。（没用过）
上代码
1 import cv2 2 import numpy as np 3 import matplotlib.pyplot as plt 4 5 6 img=cv2.imread("images/luxun.png",cv2.IMREAD_GRAYSCALE) # 读入图片 7 sobel_x = np.array([[-1, 0, 1],[-2,0,+2],[-1, 0, 1]]) # sobel的x方向算子 8 sobel_y = np.array([[1, 2, 1],[0,0,0],[-1, -2, -1]]) # sobel的x方向算子 9 sobel_x=cv2.flip(sobel_x,-1) # cv2.filter2D()计算的是相关，真正的卷积需要翻转，在进行相关计算。 10 sobel_x=cv2.flip(sobel_y,-1) 11 # cv2.flip()第二个参数：等于0：沿着x轴反转。大于0：沿着y轴反转。小于零：沿着x轴，y轴同时反转 12 13 # 卷积 opencv是用滤波器函数实现的 14 img_x=cv2.filter2D(img,-1, sobel_x) 15 img_y=cv2.filter2D(img,-1, sobel_y) 16 # 画图 plt不支持中文，但是可以通过以下方法设置修复 17 plt.rcParams['font.sans-serif']=['SimHei'] 18 plt.rcParams['axes.unicode_minus'] = False 19 20 plt.subplot(221), plt.imshow(img_x, 'gray'),plt.title('sobel_x') 21 plt.subplot(222), plt.imshow(img_y, 'gray'),plt.title('sobel_y') 22 plt.subplot(223), plt.imshow(img_y+img_x, 'gray'),plt.title('sobel') 23 plt.subplot(224), plt.imshow(img, 'gray'),plt.title('原图') 24 plt.show()
运行效果：
需要注意一点：在图像处理领域，卷积运算的定义是先将核关于x轴和y轴反转，然在做相关运算。然而工程实践中往往跳过反转，用相关运算代替卷积（比如opencv）。如果你需要严格的卷积运算，应该注意原函数的具体实现方式。sobel算子天生关于中心对称，所以反转与否并不影响结果（我在代码里用cv2.flip（）进行了反转操作）。
在之后的实现中，我发现用opencv自带的滤波函数算出来的梯度是归一化到（0-255）的，引入其他的库也很麻烦，因此自己写了个简单的二位卷积函数来实现梯度计算。所以上面的图适合看效果，并不适合在程序中使用，卷积函数的代码如下：
1 def conv2d(src,kernel): # 输入必须为方形卷积核 2 # 本函数仍然是相关运算，没有反转。如果非要严格的卷积运算，把下面一行代码的注释取消。 3 #kernel=cv2.flip(kernel,-1) 4 [rows,cols] = kernel.shape 5 border=rows//2 # 向下取整 获得卷积核边长 6 [rows,cols]=src.shape 7 dst = np.zeros(src.shape) # 采用零填充再卷积，卷积结果不会变小。 8 # print("图像长：",rows,"宽：",cols,"核边界",border) 9 # print(border,rows-border,border,cols-border) 10 temp=[] 11 for i in range(border,rows-border): 12 for j in range(border,cols-border): 13 temp=src[i-border:i+border+1,j-border:j+border+1] # 从图像获取与核匹配的图像 14 # 切片语法：索引位置包括开头但不包括结尾 [start: end: step] 15 dst[i][j]=(kernel*temp).sum() # 计算卷积 16 return dst
小技巧：用plt显示二维矩阵，鼠标移到某个像素就会显示坐标（x，y）和灰度值，浮点数也可以显示。这可以很方便的看某个数据（像素点）是否有问题。
梯度和幅值的计算效果如下：
能看出来sobel算子计算的边缘很粗很亮，比较明显，但是不够精确，我们的目标是精确到一个像素宽，至于梯度相位就很难看出什么特征，并且梯度相位实际上是为了下一步打基础的。下面附上代码：
1 img_x=conv2d(img,sobel_x) # 使用我自己的写的卷积计算梯度 2 img_y=conv2d(img,sobel_y) 3 G=np.sqrt(img_x*img_x+img_y*img_y) # 梯度幅值 4 theta=np.arctan(img_y,(img_x+0.0000000000001))*180/np.pi # 化为角度，分母+极小值是为了避免除以0 5 # plt.imshow(theta, 'gray'),plt.title('梯度相位') 6 plt.imshow(G, 'gray'),plt.title('梯度幅值') 7 plt.show()
step3:非极大值抑制
sobel算子检测出来的边缘太粗了，我们需要抑制那些梯度不够大的像素点，只保留最大的梯度，从而达到瘦边的目的。这些梯度不够大的像素点很可能是某一条边缘的过渡点。按照高数上二位函数的极大值的定义，即对点（x0，y0）的某个邻域内所有（x，y）都有f（x，y）≤(f（x0，y0），则称f在（x0，y0）具有一个极大值，极大值为f（x0，y0）。简单方案是判断一个像素点的8邻域与中心像素谁更大，但这很容易筛选出噪声，因此我们需要用梯度和梯度方向来辅助确定。
如下图所示，中心像素C的梯度方向是蓝色直线，那么只需比较中心点C与dTmp1和dTmp2的大小即可。由于这两个点的像素不知道，假设像素变化是连续的，就可以用g1、g2和g3、g4进行线性插值估计。设g1的幅值M(g1)，g2的幅值M(g2),则M(dtmp1)=w*M(g2)+(1-w)*M(g1) ，其w=distance(dtmp1,g2)/distance(g1,g2) 。也就是利用g1和g2到dTmp1的距离作为权重，来估计dTmp1的值。w在程序中可以表示为tan（θ）来表示，具体又分为四种情况（下面右图）讨论。
如下图，经过非极大值抑制可以很明显的看出去除了很多点，边缘也变得很细。在程序实现中，要注意opencv的默认坐标系是从左到右为x轴，从上到下是y轴，原点位于左上方，计算g1、g2、g3、g4的位置的时候，一定要小心（坑了我很久）。经过非极大值抑制可以看出来图片的边缘明显变细，很多看起来黑色的部分其实有值的，只是因为值太小了看不清楚，而这些黑色的部分可能是噪声或者其他原因造成的局部极大值，下一步我们就要用双阈值来限定出强边缘和弱边缘，尽可能的减少噪声的检出。代码附上：
1 # step3:非极大值抑制 2 anchor=np.where(G!=0) # 获取非零梯度的位置 3 Mask=np.zeros(img.shape) 4 5 for i in range(len(anchor[0])): 6 x=anchor[0][i] 7 y=anchor[1][i] 8 center_point=G[x,y] 9 current_theta=theta[x,y] 10 dTmp1=0 11 dTmp2=0 12 W=0 13 if current_theta>=0 and current_theta<45: 14 # g1 第一种情况 15 # g4 C g2 16 # g3 17 g1 = G[x + 1, y - 1] 18 g2 = G[x + 1, y] 19 g3 = G[x - 1, y + 1] 20 g4 = G[x - 1, y] 21 W=abs(np.tan(current_theta*np.pi/180)) # tan0-45范围为0-1 22 dTmp1= W*g1+(1-W)*g2 23 dTmp2= W*g3+(1-W)*g4 24 25 elif current_theta>=45 and current_theta<90: 26 # g2 g1 第二种情况 27 # C 28 # g3 g4 29 30 g1 = G[x + 1, y - 1] 31 g2 = G[x, y - 1] 32 g3 = G[x - 1, y + 1] 33 g4 = G[x, y + 1] 34 W = abs(np.tan((current_theta-90) * np.pi / 180)) 35 dTmp1= W*g1+(1-W)*g2 36 dTmp2= W*g3+(1-W)*g4 37 38 elif current_theta>=-90 and current_theta<-45: 39 # g1 g2 第三种情况 40 # C 41 # g4 g3 42 g1 = G[x - 1, y - 1] 43 g2 = G[x, y - 1] 44 g3 = G[x + 1, y + 1] 45 g4 = G[x, y + 1] 46 W = abs(np.tan((current_theta-90) * np.pi / 180)) 47 dTmp1= W*g1+(1-W)*g2 48 dTmp2= W*g3+(1-W)*g4 49 50 elif current_theta>=-45 and current_theta<0: 51 # g3 第四种情况 52 # g4 C g2 53 # g1 54 g1 = G[x + 1, y + 1] 55 g2 = G[x + 1, y] 56 g3 = G[x - 1, y - 1] 57 g4 = G[x - 1, y] 58 W = abs(np.tan(current_theta * np.pi / 180)) 59 dTmp1= W*g1+(1-W)*g2 60 dTmp2= W*g3+(1-W)*g4 61 62 if dTmp1<center_point and dTmp2<center_point: # 记录极大值结果 63 Mask[x,y]=center_point 64 #Mask=(Mask-Mask.min())/(Mask.max()-Mask.min())*256 #归一化 65 plt.imshow(Mask,'gray'),plt.title('Mask') 66 plt.show()
step4：用双阈值算法检测和连接边缘
双阈值法非常简单，我们假设两类边缘：经过非极大值抑制之后的边缘点中，梯度值超过T1的称为强边缘，梯度值小于T1大于T2的称为弱边缘，梯度小于T2的不是边缘。可以肯定的是，强边缘必然是边缘点，因此必须将T1设置的足够高，以要求像素点的梯度值足够大（变化足够剧烈），而弱边缘可能是边缘，也可能是噪声，如何判断呢？当弱边缘的周围8邻域有强边缘点存在时，就将该弱边缘点变成强边缘点，以此来实现对强边缘的补充。实际中人们发现T1:T2=2:1的比例效果比较好，其中T1可以人为指定，也可以设计算法来自适应的指定，比如定义梯度直方图的前30%的分界线为T1，我实现的是人为指定阈值。检查8邻域的方法叫边缘滞后跟踪，连接边缘的办法还有区域生长法等等。强边缘、弱边缘、综合效果、和opencv的canny函数对比如下：
三、总结
实现结果还是很打击的，我检测到的边缘过于断续，没有opencv实现的效果好。查了一下opencv的源码，这里猜测两个可能的原因：源码里梯度的方向被近似到四个角度之一 （0，45，90，135），但我用线性插值的的结果是梯度方向更精确，而过于精确-->过于严格-->容易受到噪声干扰，所以在非极大值抑制这之后，我比opencv少了更多的点，最终导致了边缘不够连续；第二个原因可能是边缘连接算法效果不够好，把图象放大来看，我产生的边缘倾向于对角线上连接，而opencv的边缘倾向于折线连接，因此opencv的边缘更完整连续，而我的边缘更细，更容易断续。
限于时间，暂时研究到这里，希望各位多多指正！感谢所有我参考过的博客和文档！
1 import cv2 2 import numpy as np 3 import matplotlib.pyplot as plt 4 5 # 画图 plt不支持中文，但是可以通过以下方法设置修复 6 plt.rcParams['font.sans-serif']=['SimHei'] 7 plt.rcParams['axes.unicode_minus'] = False 8 9 def conv2d(src,kernel): # 输入必须为方形卷积核 10 # 本函数仍然是相关运算，没有反转。如果非要严格的卷积运算，把下面一行代码的注释取消。 11 #kernel=cv2.flip(kernel,-1) 12 [rows,cols] = kernel.shape 13 border=rows//2 # 向下取整 获得卷积核边长 14 [rows,cols]=src.shape 15 dst = np.zeros(src.shape) # 采用零填充再卷积，卷积结果不会变小。 16 # print("图像长：",rows,"宽：",cols,"核边界",border) 17 # print(border,rows-border,border,cols-border) 18 temp=[] 19 for i in range(border,rows-border): 20 for j in range(border,cols-border): 21 temp=src[i-border:i+border+1,j-border:j+border+1] # 从图像获取与核匹配的图像 22 # 切片语法：索引位置包括开头但不包括结尾 [start: end: step] 23 dst[i][j]=(kernel*temp).sum() # 计算卷积 24 return dst 25 26 27 # step0:读入图片 28 img=cv2.imread("images/luxun.png",cv2.IMREAD_GRAYSCALE) # 读入图片 29 30 # step1:高斯滤波 31 img=cv2.GaussianBlur(img,(5,5),0) 32 33 # step2:计算梯度强度和方向 34 sobel_x = np.array([[-1, 0, 1],[-2,0,+2],[-1, 0, 1]]) # sobel的x方向算子 35 sobel_y = np.array([[1, 2, 1],[0,0,0],[-1, -2, -1]]) # sobel的y方向算子 36 37 # img_x=cv2.filter2D(img,-1, sobel_x) # 这个滤波器会将卷积结果归一化到0-255，无法计算梯度方向。 38 # img_y=cv2.filter2D(img,-1, sobel_y) # 而真正的图像卷积可能会出现负数，因此只能自己写个卷积。 39 img_x=conv2d(img,sobel_x) # 使用我自己的写的卷积计算梯度 40 img_y=conv2d(img,sobel_y) 41 G=np.sqrt(img_x*img_x+img_y*img_y) # 梯度幅值 42 theta=np.arctan(img_y,(img_x+0.0000000000001))*180/np.pi # 化为角度，分母+极小值是为了避免除以0 43 44 # plt.imshow(theta, 'gray'),plt.title('梯度相位') 45 # plt.imshow(G, 'gray'),plt.title('梯度幅值') 46 # plt.show() 47 # exit() 48 49 # step3:非极大值抑制 50 anchor=np.where(G!=0) # 获取非零梯度的位置 51 Mask=np.zeros(img.shape) 52 for i in range(len(anchor[0])): 53 x=anchor[0][i] # 取出第i个非零梯度的x坐标 54 y=anchor[1][i] 55 center_point=G[x,y] 56 current_theta=theta[x,y] 57 dTmp1=0 58 dTmp2=0 59 W=0 60 if current_theta>=0 and current_theta<45: 61 # g1 第一种情况 62 # g4 C g2 63 # g3 64 g1 = G[x + 1, y - 1] 65 g2 = G[x + 1, y] 66 g3 = G[x - 1, y + 1] 67 g4 = G[x - 1, y] 68 W=abs(np.tan(current_theta*np.pi/180)) # tan0-45范围为0-1 69 dTmp1= W*g1+(1-W)*g2 70 dTmp2= W*g3+(1-W)*g4 71 72 elif current_theta>=45 and current_theta<90: 73 # g2 g1 第二种情况 74 # C 75 # g3 g4 76 g1 = G[x + 1, y - 1] 77 g2 = G[x, y - 1] 78 g3 = G[x - 1, y + 1] 79 g4 = G[x, y + 1] 80 W = abs(np.tan((current_theta-90) * np.pi / 180)) 81 dTmp1= W*g1+(1-W)*g2 82 dTmp2= W*g3+(1-W)*g4 83 elif current_theta>=-90 and current_theta<-45: 84 # g1 g2 第三种情况 85 # C 86 # g4 g3 87 g1 = G[x - 1, y - 1] 88 g2 = G[x, y - 1] 89 g3 = G[x + 1, y + 1] 90 g4 = G[x, y + 1] 91 W = abs(np.tan((current_theta-90) * np.pi / 180)) 92 dTmp1= W*g1+(1-W)*g2 93 dTmp2= W*g3+(1-W)*g4 94 elif current_theta>=-45 and current_theta<0: 95 # g3 第四种情况 96 # g4 C g2 97 # g1 98 g1 = G[x + 1, y + 1] 99 g2 = G[x + 1, y] 100 g3 = G[x - 1, y - 1] 101 g4 = G[x - 1, y] 102 W = abs(np.tan(current_theta * np.pi / 180)) 103 dTmp1= W*g1+(1-W)*g2 104 dTmp2= W*g3+(1-W)*g4 105 106 if dTmp1<center_point and dTmp2<center_point: # 记录极大值结果 107 Mask[x,y]=center_point 108 # plt.imshow(Mask,'gray'),plt.title('Mask') 109 # plt.show() 110 # exit() 111 112 # step4:双阈值选取 113 high_threshold=100 114 low_threshold=high_threshold/2 115 strong_edge=np.zeros(G.shape) # 强边缘 116 weak_edge=np.zeros(G.shape) # 弱边缘 117 118 xNum = [1, 1, 0, -1, -1, -1, 0, 1] # 8邻域偏移坐标 119 yNum = [0, 1, 1, 1, 0, -1, -1, -1] 120 [rows, cols] = G.shape 121 for i in range(rows): 122 for j in range(cols): 123 current_point=Mask[i,j] 124 if current_point>0: 125 if current_point>high_threshold: # 强边缘提取 126 strong_edge[i,j]=255 127 elif current_point<high_threshold and current_point>low_threshold: # 弱边缘提取 128 # step6:顺便进行边缘连接 129 change = True 130 while change: 131 change = False 132 for k in range(8): 133 xx=i+xNum[k] 134 yy=j+yNum[k] 135 if Mask[xx,yy]>high_threshold: 136 weak_edge[i, j] = 255 137 break # 跳出八邻域循环 138 output=strong_edge+weak_edge # 强弱边缘综合效果 139 140 img_edge = cv2.Canny(img, 50, 100) # opencv实现效果 141 142 # 显示效果 143 plt.subplot(141), plt.imshow(strong_edge, 'gray'),plt.title('strong_edge') 144 plt.subplot(142), plt.imshow(weak_edge, 'gray'),plt.title('weak_edge') 145 plt.subplot(143), plt.imshow(output, 'gray'),plt.title('result') 146 plt.subplot(144), plt.imshow(img_edge, 'gray'),plt.title('opencv') 147 plt.show()
参考文献：
https://www.cnblogs.com/love6tao/p/5152020.html
https://www.cnblogs.com/techyan1990/p/7291771.html