计算机视觉学习笔记#2 边缘提取:Canny(附代码)

本文最后更新于:2022年7月21日 下午

计算机视觉学习笔记#2 边缘提取(附代码)

边缘直观来说指的是图像中突然的或者非连续的变化,提取到了图像中的边缘后能有助于理解图像。如下图,了解了边缘后就能明白图片是一个女人。

理想的边缘

边缘的产生

现实当中,边缘的产生可能因为:

  • 面的不连续
  • 深度的不连续
  • 表面颜色的不连续
  • 光照带来的不连续

边缘的产生

用导数获取到边缘

这一段在数字图像处理中已经学过所以只是简单描述:

一条有宽度的线如图,根据红线可以得到一条强度函数,然后求一阶导数,发现只需要关注极值点就能找到边缘。

边缘的数学表示

但是直接用[-1,1]或者[-1;1]这样的卷积核太简单了,所以之后出现了Prewitt、Sobel、Robert算子。同时,使用$||\nabla f||=\sqrt{M_x^2+M_y^2}$可以求出强度,用$\theta=tan^{-1}(M_x/M_y)$可以求梯度方向。

若要考试,要注意0的方向是检测线的方向,-1 0 1的方向和线的方向垂直。

各类边缘检测算子

在实际情况中,可能有非常多的噪声干扰边缘提取,会出现很多伪边缘,所以一般先进行一个高斯平滑再求导:$\frac{d}{dx}(f*g)$,其中$g$是平滑模板。

而上式等于$f*\frac{d}{dx}g$,所以只需要卷积一次$\frac{d}{dx}g$就可以了。

此时,高斯平滑带来了一个超参数方差$\sigma$和窗大小,而窗经常设置成$\pm3\sigma$。

高斯卷积模版

高斯核值之和为1,不带来信号放大;而高斯梯度核之和为0,平坦地区无响应。高斯核没有负数,高斯梯度核有负数。

对比高斯核和高斯梯度核

Canny边缘提取方法

Canny是一种有名的边缘提取方法,基于我们上面找到的方法,Canny先用高斯梯度核卷积一遍,然后求$||\nabla f||=\sqrt{M_x^2+M_y^2}$。

之后会进行非最大化抑制:即对于每一个像素,在梯度方向上观察周围的像素中自己是否为最大值,假如是则保留,假如不是则置0。

最后使用双门限过滤,先用高门限得到更少但是更强的边缘,然后用低门限得到更多但是噪声也多的边缘,然后认为低门限中能连接高门限两条边的边缘也是有用的边缘,添加回去。

Canny

一份手写的Canny算法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
import numpy as np
import cv2
import matplotlib.pyplot as plt


def get_gaussian(size=5, sigma=1):
assert sigma > 0 and size >= 3, "参数错误"
x, y = np.meshgrid(
np.linspace(-3 * sigma, 3 * sigma, size),
np.linspace(-3 * sigma, 3 * sigma, size)
)
gaussian = 1 / (2 * np.pi * sigma ** 2) * np.exp(-(x ** 2 + y ** 2) / (2 * sigma ** 2))
return gaussian / gaussian.sum()


def simple_2d_conv(x: np.ndarray, kernel: np.ndarray, mode: str = 'same'):
"""
对灰度图像进行2d卷积操作
:param x: 输入图像
:param kernel: 卷积核
:param mode: 卷积模式 full same valid
:return:
"""
assert len(x.shape) == 2 and len(kernel.shape) == 2, "只支持灰度图像"
assert kernel.shape[0] % 2 == 1 and kernel.shape[1] % 2 == 1, "参数错误"
assert mode in ['same', 'full', 'valid'], "卷积模式错误"
h, w = x.shape # h行 w列
kh, kw = kernel.shape
# padding
if mode == 'full':
# np.pad第二个参数是 上下左右 的顺序
x = np.pad(x, ((kh - 1, kh - 1), (kw - 1, kw - 1)), 'constant', constant_values=0)
h += kh - 1
w += kw - 1
elif mode == 'same':
# np.pad第二个参数是 上下左右 的顺序
x = np.pad(x, ((kh // 2, kh // 2), (kw // 2, kw // 2)), 'constant', constant_values=0)
else:
h -= kh - 1
w -= kw - 1
# conv
res = np.zeros([h, w])
for i in range(h):
for j in range(w):
res[i, j] = np.sum(x[i:i + kh, j:j + kw] * kernel)
return res


def non_max_suppression(x, theta):
res = np.zeros_like(x)
for i in range(1, x.shape[0] - 1):
for j in range(1, x.shape[1] - 1):
abs_theta = abs(theta[i, j])
# 分成4个区域求插值
if abs_theta <= np.pi / 4:
interpolation1 = np.tan(abs_theta) * x[i + 1, j] + (1 - np.tan(abs_theta)) * x[i + 1, j + 1]
interpolation2 = np.tan(abs_theta) * x[i - 1, j] + (1 - np.tan(abs_theta)) * x[i - 1, j - 1]
elif abs_theta <= np.pi / 2:
interpolation1 = 1 / np.tan(abs_theta) * x[i, j + 1] + (1 - 1 / np.tan(abs_theta)) * x[i + 1, j + 1]
interpolation2 = 1 / np.tan(abs_theta) * x[i, j - 1] + (1 - 1 / np.tan(abs_theta)) * x[i - 1, j - 1]
elif abs_theta <= np.pi * 3 / 4:
interpolation1 = np.tan(abs_theta) * x[i, j + 1] + (1 - np.tan(abs_theta)) * x[i - 1, j + 1]
interpolation2 = np.tan(abs_theta) * x[i, j - 1] + (1 - np.tan(abs_theta)) * x[i + 1, j - 1]
else:
interpolation1 = -np.tan(abs_theta) * x[i - 1, j] + (1 + np.tan(abs_theta)) * x[i - 1, j + 1]
interpolation2 = -np.tan(abs_theta) * x[i + 1, j] + (1 + np.tan(abs_theta)) * x[i + 1, j - 1]
if x[i, j] >= max(interpolation1, interpolation2):
res[i, j] = x[i, j]
return res


def double_thresholding(x, low_threshold, high_threshold):
# 1st-pass 标记所有点
thres_map = np.zeros_like(x)
for i in range(x.shape[0]):
for j in range(x.shape[1]):
if x[i, j] >= high_threshold:
thres_map[i, j] = 2
elif x[i, j] >= low_threshold:
thres_map[i, j] = 1
# 2nd-pass 弱边界进行连接的判断
# 和OpenCV的实现差距可能就在这里
cc = 1
for i in range(cc, x.shape[0] - cc):
for j in range(cc, x.shape[1] - cc):
if thres_map[i, j] == 1:
if np.any(thres_map[i - cc:i + cc, j - cc:j + cc] == 2):
thres_map[i, j] = 2
else:
thres_map[i, j] = 0
thres_map[thres_map == 1] = 0
return thres_map


img = cv2.imread("tiger.jpg", cv2.IMREAD_GRAYSCALE)
# 1. 高斯平滑
# kernel = get_gaussian(size=5, sigma=1)
kernel = np.array([[2, 4, 5, 4, 2], [4, 9, 12, 9, 4], [5, 12, 15, 12, 5], [4, 9, 12, 9, 4], [2, 4, 5, 4, 2]]) / 159
conv_img = simple_2d_conv(img, kernel, 'valid')
# 2. Sobel滤波
Sobel_X = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])
Sobel_Y = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
sobel_x_img = simple_2d_conv(conv_img, Sobel_X, 'valid')
sobel_y_img = simple_2d_conv(conv_img, Sobel_Y, 'valid')
intensity = np.abs(sobel_x_img) + np.abs(sobel_y_img)
intensity = intensity / np.max(intensity) * 255
# 3. 非极大抑制
theta_map = np.arctan2(sobel_y_img, sobel_x_img) # [-pi, +pi]
non_max_img = non_max_suppression(intensity, theta_map)
# 4. 双阈值
# print(np.max(non_max_img), np.min(non_max_img))
double_threshold_img = double_thresholding(non_max_img, 5, 40)

# show
plt.subplot(2, 2, 1)
plt.imshow(img, cmap=plt.get_cmap('gray'))
plt.title('Original')
plt.subplot(2, 2, 2)
plt.imshow(conv_img, cmap=plt.get_cmap('gray'))
plt.title('Gaussian Filter')
plt.subplot(2, 2, 3)
plt.imshow(cv2.Canny(img, 90, 200), cmap=plt.get_cmap('gray'))
plt.title('OpenCV Canny')
plt.subplot(2, 2, 4)
plt.imshow(double_threshold_img, cmap=plt.get_cmap('gray'))
plt.title('Canny')

plt.show()