最优边缘准则

Canny 的目标是找到一个最优的边缘检测算法，最优边缘检测的含义是：
(1)最优检测：算法能够尽可能多地标识出图像中的实际边缘，漏检真实边缘的概率和误检非边缘的概率都尽可能小；
(2)最优定位准则：检测到的边缘点的位置距离实际边缘点的位置最近，或者是由于噪声影响引起检测出的边缘偏离物体的真实边缘的程度最小；
(3)检测点与边缘点一一对应：算子检测的边缘点与实际边缘点应该是一一对应。

算法实现步骤

Canny边缘检测算法可以分为以下5个步骤：

1. 应用高斯滤波来平滑(模糊)图像，目的是去除噪声

高斯滤波器是将高斯函数离散化，将滤波器中对应的横纵坐标索引代入到高斯函数，从而得到对应的值。

二维的高斯函数如下：其中 (x , y)为坐标， σ 为标准差
$\frac{1}{2\pi σ^2} e^{- \frac{x^2 + y^2}{2σ^2}} \tag1$
不同尺寸的滤波器，得到的值也不同，下面是 (2k+1)x(2k+1) 滤波器的计算公式 :
$\frac{1}{2\pi σ^2} e^{- \frac{(i-k-1)^2 + (j-k-1)^2}{2σ^2}} \tag2$
常见的高斯滤波器大小为 5*5， σ = 1.4 ，其近似值为：

2. 计算梯度强度和方向

接下来，我们要寻找边缘，即灰度强度变化最强的位置，（一道黑边一道白边中间就是边缘，它的灰度值变化是最大的）。在图像中，用梯度来表示灰度值的变化程度和方向。
常见方法采用Sobel滤波器【水平x和垂直y方向】在计算梯度和方向
水平方向的Sobel算子G_x：用来检测 y 方向的边缘

-1	0	1
-2	0	2
-1	0	1

垂直方向的Sobel算子G_y：用来检测 x 方向的边缘（ 边缘方向和梯度方向垂直）

1	2	1
0	0	0
-1	-2	-1

采用下列公式计算梯度和方向：
$\sqrt{(G_x^2 + G_y^2)} \tag3$
$\theta = arctan{\frac{G_y}{G_x}} \tag4$

3. 应用非最大抑制技术NMS来消除边误检

原理：遍历梯度矩阵上的所有点，并保留边缘方向上具有极大值的像素

这一步的目的是将模糊（blurred）的边界变得清晰（sharp）。通俗的讲，就是保留了每个像素点上梯度强度的极大值，而删掉其他的值。对于每个像素点，进行如下操作：
a) 将其梯度方向近似为以下值中的一个（0,45,90,135,180,225,270,315）（即上下左右和45度方向）

b) 比较该像素点，和其梯度方向正负方向的像素点的梯度强度

c) 如果该像素点梯度强度最大则保留，否则抑制（删除，即置为0）
$M_T(m,n) = \begin{cases} M(m,n), & \text {if M(m,n) > T}\\ 0, & \text {otherwise} \end{cases}$
例如：【该例子来自 Python - Opencv 之 Canny 边缘检测】

点 A 位于图像边缘垂直方向. 梯度方向垂直于边缘. 点 B 和点 C 位于梯度方向. 因此，检查点 A 和点 B，点 C，确定点A是否是局部最大值. 如果点 A 是局部最大值，则继续下一个阶段；如果点 A 不是局部最大值，则其被抑制设为0。

最后会保留一条边界处最亮的一条细线

4. 应用双阈值的方法来决定可能的（潜在的）边界

这个阶段决定哪些边缘是真正的边缘，哪些边缘不是真正的边缘

经过非极大抑制后图像中仍然有很多噪声点。Canny算法中应用了一种叫双阈值的技术。即设定一个阈值上界maxVal和阈值下界minVal，图像中的像素点如果大于阈值上界则认为必然是边界（称为强边界，strong edge），小于阈值下界则认为必然不是边界，两者之间的则认为是候选项（称为弱边界，weak edge），需进行进一步处理——如果与确定为边缘的像素点邻接，则判定为边缘；否则为非边缘。

5. 利用滞后技术来跟踪边界

这个阶段是进一步处理弱边界

大体思想是，和强边界相连的弱边界认为是边界，其他的弱边界则被抑制。
由真实边缘引起的弱边缘像素将连接到强边缘像素，而噪声响应未连接。为了跟踪边缘连接，通过查看弱边缘像素及其8个邻域像素，只要其中一个为强边缘像素，则该弱边缘点就可以保留为真实的边缘。

opencv实现Canny边缘检测

OpenCV 提供了 cv2.canny 函数.

edge = cv2.Canny(image, threshold1, threshold2[, edges[, apertureSize[, L2gradient ]]])

参数 image - 输入图片，必须为单通道的灰度图
参数 threshold1 和 threshold2 - 分别对应于阈值 minVal 和 maxVal
参数 apertureSize - 用于计算图片提取的 Sobel kernel 尺寸. 默认为 3.
参数 L2gradient - 指定计算梯度的等式. 当参数为 True 时，采用梯度计算公式(3)(4)，其精度更高；否则采用的梯度计算公式为： $G = |G_x| + |G_y|$ . 该参数默认为 False.

e.g.

import numpy as np
import cv2 as cv
from matplotlib import pyplot as pltimg = cv.imread('test.jpg',0)
edges = cv.Canny(img, 100, 200)plt.subplot(121),plt.imshow(img,cmap = 'gray')
plt.title('Original Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(edges,cmap = 'gray')
plt.title('Edge Image'), plt.xticks([]), plt.yticks([])plt.show()

结果如下图：
Canny边缘检测结果

手写代码

声明：此部分代码【from Python实现Canny算子边缘检测】

import numpy as np
import cv2 as cv
from matplotlib import pyplot as pltdef smooth(image, sigma = 1.4, length = 5):""" Smooth the imageCompute a gaussian filter with sigma = sigma and kernal_length = length.Each element in the kernal can be computed as below:G[i, j] = (1/(2*pi*sigma**2))*exp(-((i-k-1)**2 + (j-k-1)**2)/2*sigma**2)Then, use the gaussian filter to smooth the input image.Args:image: array of grey imagesigma: the sigma of gaussian filter, default to be 1.4length: the kernal length, default to be 5Returns:the smoothed image"""# Compute gaussian filterk = length // 2gaussian = np.zeros([length, length])for i in range(length):for j in range(length):gaussian[i, j] = np.exp(-((i-k) ** 2 + (j-k) ** 2) / (2 * sigma ** 2))gaussian /= 2 * np.pi * sigma ** 2# Batch Normalizationgaussian = gaussian / np.sum(gaussian)# Use Gaussian FilterW, H = image.shapenew_image = np.zeros([W - k * 2, H - k * 2])for i in range(W - 2 * k):for j in range(H - 2 * k):# 卷积运算new_image[i, j] = np.sum(image[i:i+length, j:j+length] * gaussian)new_image = np.uint8(new_image)return new_imagedef get_gradient_and_direction(image):""" Compute gradients and its directionUse Sobel filter to compute gradients and direction.-1 0 1        -1 -2 -1Gx = -2 0 2   Gy =  0  0  0-1 0 1         1  2  1Args:image: array of grey imageReturns:gradients: the gradients of each pixeldirection: the direction of the gradients of each pixel"""Gx = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])Gy = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])W, H = image.shapegradients = np.zeros([W - 2, H - 2])direction = np.zeros([W - 2, H - 2])for i in range(W - 2):for j in range(H - 2):dx = np.sum(image[i:i+3, j:j+3] * Gx)dy = np.sum(image[i:i+3, j:j+3] * Gy)gradients[i, j] = np.sqrt(dx ** 2 + dy ** 2)if dx == 0:direction[i, j] = np.pi / 2else:direction[i, j] = np.arctan(dy / dx)gradients = np.uint8(gradients)return gradients, directiondef NMS(gradients, direction):""" Non-maxima suppressionArgs:gradients: the gradients of each pixeldirection: the direction of the gradients of each pixelReturns:the output image"""W, H = gradients.shapenms = np.copy(gradients[1:-1, 1:-1])for i in range(1, W - 1):for j in range(1, H - 1):theta = direction[i, j]weight = np.tan(theta)if theta > np.pi / 4:d1 = [0, 1]d2 = [1, 1]weight = 1 / weightelif theta >= 0:d1 = [1, 0]d2 = [1, 1]elif theta >= - np.pi / 4:d1 = [1, 0]d2 = [1, -1]weight *= -1else:d1 = [0, -1]d2 = [1, -1]weight = -1 / weightg1 = gradients[i + d1[0], j + d1[1]]g2 = gradients[i + d2[0], j + d2[1]]g3 = gradients[i - d1[0], j - d1[1]]g4 = gradients[i - d2[0], j - d2[1]]grade_count1 = g1 * weight + g2 * (1 - weight)grade_count2 = g3 * weight + g4 * (1 - weight)if grade_count1 > gradients[i, j] or grade_count2 > gradients[i, j]:nms[i - 1, j - 1] = 0return nmsdef double_threshold(nms, threshold1, threshold2):""" Double ThresholdUse two thresholds to compute the edge.Args:nms: the input imagethreshold1: the low thresholdthreshold2: the high thresholdReturns:The binary image."""visited = np.zeros_like(nms)output_image = nms.copy()W, H = output_image.shapedef dfs(i, j):if i >= W or i < 0 or j >= H or j < 0 or visited[i, j] == 1:returnvisited[i, j] = 1if output_image[i, j] > threshold1:output_image[i, j] = 255dfs(i-1, j-1)dfs(i-1, j)dfs(i-1, j+1)dfs(i, j-1)dfs(i, j+1)dfs(i+1, j-1)dfs(i+1, j)dfs(i+1, j+1)else:output_image[i, j] = 0for w in range(W):for h in range(H):if visited[w, h] == 1:continueif output_image[w, h] >= threshold2:dfs(w, h)elif output_image[w, h] <= threshold1:output_image[w, h] = 0visited[w, h] = 1for w in range(W):for h in range(H):if visited[w, h] == 0:output_image[w, h] = 0return output_imageif __name__ == "__main__":# code to read imageimage = cv.imread('test.jpg',0)cv.imshow("Original",image)smoothed_image = smooth(image)cv.imshow("GaussinSmooth(5*5)",smoothed_image)gradients, direction = get_gradient_and_direction(smoothed_image)# print(gradients)# print(direction)nms = NMS(gradients, direction)output_image = double_threshold(nms, 40, 100)cv.imshow("outputImage",output_image)cv.waitKey(0)