实现Canny算法原理 python

Posted 2021-03-04 我坚信阳光灿烂

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了实现Canny算法原理 python相关的知识，希望对你有一定的参考价值。

创作很累，如果您觉得对您有帮助，请点赞支持，感谢！

一. 总的算法流程：

① 使用高斯滤波器滤波

② 使用 Sobel 滤波器滤波获得在 x 和 y 方向上的输出，在此基础上求出边缘的强度和边缘的角度

edge 为边缘强度，tan 为边缘角度 ↑

③ 对边缘角度进行量化处理

对边缘角度进行量化处理算法 ↑

④ 根据边缘角度对边缘强度进行非极大值抑制（Non-maximum suppression），使图像边缘变得更细

非极大值抑制算法：0°时取(x,y)、(x+1,y)、(x-1,y) 中的最大值，其它角度类似 ↑

⑤ 使用滞后阈值对图像进行二值化处理，优化图像显示效果

算法如上所示 ↑

⑥ 输出图像边缘提取效果

二. 使用python手动实现 Canny 算法，完成图像边缘提取

  1 # writer：wojianxinygcl@163.com
  2 # date  ：2020.3.20
  3 import cv2
  4 import numpy as np
  5 import matplotlib.pyplot as plt
  6 
  7 def Canny(img):
  8 
  9     # Gray scale
 10     def BGR2GRAY(img):
 11         b = img[:, :, 0].copy()
 12         g = img[:, :, 1].copy()
 13         r = img[:, :, 2].copy()
 14 
 15         # Gray scale
 16         out = 0.2126 * r + 0.7152 * g + 0.0722 * b
 17         out = out.astype(np.uint8)
 18 
 19         return out
 20 
 21 
 22     # Gaussian filter for grayscale
 23     def gaussian_filter(img, K_size=3, sigma=1.4):
 24 
 25         if len(img.shape) == 3:
 26             H, W, C = img.shape
 27             gray = False
 28         else:
 29             img = np.expand_dims(img, axis=-1)
 30             H, W, C = img.shape
 31             gray = True
 32 
 33         ## Zero padding
 34         pad = K_size // 2
 35         out = np.zeros([H + pad * 2, W + pad * 2, C], dtype=np.float)
 36         out[pad : pad + H, pad : pad + W] = img.copy().astype(np.float)
 37 
 38         ## prepare Kernel
 39         K = np.zeros((K_size, K_size), dtype=np.float)
 40         for x in range(-pad, -pad + K_size):
 41             for y in range(-pad, -pad + K_size):
 42                 K[y + pad, x + pad] = np.exp( - (x ** 2 + y ** 2) / (2 * sigma * sigma))
 43         #K /= (sigma * np.sqrt(2 * np.pi))
 44         K /= (2 * np.pi * sigma * sigma)
 45         K /= K.sum()
 46 
 47         tmp = out.copy()
 48 
 49         # filtering
 50         for y in range(H):
 51             for x in range(W):
 52                 for c in range(C):
 53                     out[pad + y, pad + x, c] = np.sum(K * tmp[y : y + K_size, x : x + K_size, c])
 54 
 55         out = np.clip(out, 0, 255)
 56         out = out[pad : pad + H, pad : pad + W]
 57         out = out.astype(np.uint8)
 58 
 59         if gray:
 60             out = out[..., 0]
 61 
 62         return out
 63 
 64 
 65     # sobel filter
 66     def sobel_filter(img, K_size=3):
 67         if len(img.shape) == 3:
 68             H, W, C = img.shape
 69         else:
 70             H, W = img.shape
 71 
 72         # Zero padding
 73         pad = K_size // 2
 74         out = np.zeros((H + pad * 2, W + pad * 2), dtype=np.float)
 75         out[pad : pad + H, pad : pad + W] = img.copy().astype(np.float)
 76         tmp = out.copy()
 77 
 78         out_v = out.copy()
 79         out_h = out.copy()
 80 
 81         ## Sobel vertical
 82         Kv = [[1., 2., 1.],[0., 0., 0.], [-1., -2., -1.]]
 83         ## Sobel horizontal
 84         Kh = [[1., 0., -1.],[2., 0., -2.],[1., 0., -1.]]
 85 
 86         # filtering
 87         for y in range(H):
 88             for x in range(W):
 89                 out_v[pad + y, pad + x] = np.sum(Kv * (tmp[y : y + K_size, x : x + K_size]))
 90                 out_h[pad + y, pad + x] = np.sum(Kh * (tmp[y : y + K_size, x : x + K_size]))
 91 
 92         out_v = np.clip(out_v, 0, 255)
 93         out_h = np.clip(out_h, 0, 255)
 94 
 95         out_v = out_v[pad : pad + H, pad : pad + W]
 96         out_v = out_v.astype(np.uint8)
 97         out_h = out_h[pad : pad + H, pad : pad + W]
 98         out_h = out_h.astype(np.uint8)
 99 
100         return out_v, out_h
101 
102 
103     # get edge strength and edge angle
104     def get_edge_angle(fx, fy):
105         # get edge strength
106         edge = np.sqrt(np.power(fx.astype(np.float32), 2) + np.power(fy.astype(np.float32), 2))
107         edge = np.clip(edge, 0, 255)
108 
109         # make sure the denominator is not 0
110         fx = np.maximum(fx, 1e-10)
111         #fx[np.abs(fx) <= 1e-5] = 1e-5
112 
113         # get edge angle
114         angle = np.arctan(fy / fx)
115 
116         return edge, angle
117 
118     
119     # 将角度量化为0°、45°、90°、135°
120     def angle_quantization(angle):
121         angle = angle / np.pi * 180
122         angle[angle < -22.5] = 180 + angle[angle < -22.5]
123         _angle = np.zeros_like(angle, dtype=np.uint8)
124         _angle[np.where(angle <= 22.5)] = 0
125         _angle[np.where((angle > 22.5) & (angle <= 67.5))] = 45
126         _angle[np.where((angle > 67.5) & (angle <= 112.5))] = 90
127         _angle[np.where((angle > 112.5) & (angle <= 157.5))] = 135
128 
129         return _angle
130 
131 
132     def non_maximum_suppression(angle, edge):
133         H, W = angle.shape
134         _edge = edge.copy()
135         
136         for y in range(H):
137             for x in range(W):
138                     if angle[y, x] == 0:
139                             dx1, dy1, dx2, dy2 = -1, 0, 1, 0
140                     elif angle[y, x] == 45:
141                             dx1, dy1, dx2, dy2 = -1, 1, 1, -1
142                     elif angle[y, x] == 90:
143                             dx1, dy1, dx2, dy2 = 0, -1, 0, 1
144                     elif angle[y, x] == 135:
145                             dx1, dy1, dx2, dy2 = -1, -1, 1, 1
146                     # 边界处理
147                     if x == 0:
148                             dx1 = max(dx1, 0)
149                             dx2 = max(dx2, 0)
150                     if x == W-1:
151                             dx1 = min(dx1, 0)
152                             dx2 = min(dx2, 0)
153                     if y == 0:
154                             dy1 = max(dy1, 0)
155                             dy2 = max(dy2, 0)
156                     if y == H-1:
157                             dy1 = min(dy1, 0)
158                             dy2 = min(dy2, 0)
159                     # 如果不是最大值，则将这个位置像素值置为0
160                     if max(max(edge[y, x], edge[y + dy1, x + dx1]), edge[y + dy2, x + dx2]) != edge[y, x]:
161                             _edge[y, x] = 0
162 
163         return _edge
164 
165 
166     # 滞后阈值处理二值化图像
167     # > HT 的设为255，< LT 的设置0，介于它们两个中间的值，使用8邻域判断法
168     def hysterisis(edge, HT=100, LT=30):
169         H, W = edge.shape
170 
171         # Histeresis threshold
172         edge[edge >= HT] = 255
173         edge[edge <= LT] = 0
174 
175         _edge = np.zeros((H + 2, W + 2), dtype=np.float32)
176         _edge[1 : H + 1, 1 : W + 1] = edge
177 
178         ## 8 - Nearest neighbor
179         nn = np.array(((1., 1., 1.), (1., 0., 1.), (1., 1., 1.)), dtype=np.float32)
180 
181         for y in range(1, H+2):
182                 for x in range(1, W+2):
183                         if _edge[y, x] < LT or _edge[y, x] > HT:
184                                 continue
185                         if np.max(_edge[y-1:y+2, x-1:x+2] * nn) >= HT:
186                                 _edge[y, x] = 255
187                         else:
188                                 _edge[y, x] = 0
189 
190         edge = _edge[1:H+1, 1:W+1]
191                                 
192         return edge
193 
194     # grayscale
195     gray = BGR2GRAY(img)
196 
197     # gaussian filtering
198     gaussian = gaussian_filter(gray, K_size=5, sigma=1.4)
199 
200     # sobel filtering
201     fy, fx = sobel_filter(gaussian, K_size=3)
202 
203     # get edge strength, angle
204     edge, angle = get_edge_angle(fx, fy)
205 
206     # angle quantization
207     angle = angle_quantization(angle)
208 
209     # non maximum suppression
210     edge = non_maximum_suppression(angle, edge)
211 
212     # hysterisis threshold
213     out = hysterisis(edge, 80, 20)
214 
215     return out
216 
217 
218 if __name__ == ‘__main__‘:
219     # Read image
220     img = cv2.imread("../paojie.jpg").astype(np.float32)
221 
222     # Canny
223     edge = Canny(img)
224 
225     out = edge.astype(np.uint8)
226 
227     # Save result
228     cv2.imwrite("out.jpg", out)
229     cv2.imshow("result", out)
230     cv2.waitKey(0)
231     cv2.destroyAllWindows()