Python利用Canny算法检测硬币边缘

一、问题背景

纸面上有一枚一元钱的银币，你能在 Canny 和 Hough 的帮助下找到它的坐标方程吗？

确定一个圆的坐标方程，首先我们要检测到其边缘，然后求出其在纸面上的相对位置以及半径大小。

在这篇文章中我们使用 Canny 算法来检测出纸面上银币的边缘。

二、Canny 算法

Canny 可以用于拿到图像中物体的边缘，其步骤如下

• 进行高斯平滑

• 计算图像梯度(记录其强度、方向)

• 进行非极大化抑制

• 进行滞后边缘跟踪

进行上面的四步之后，我们拿到的纸面上硬币边缘提取效果图如下

（一）、高斯平滑

class GaussianSmoothingNet(nn.Module):
   def __init__(self) -> None:
       super(GaussianSmoothingNet, self).__init__()

filter_size = 5
       # shape为(1, 5), 方差为 1.0 的高斯滤波核
       generated_filters = gaussian(filter_size,std=1.0).reshape([1,filter_size])

# GFH(V): gaussian filter of horizontal(vertical) 水平（竖直）方向的高斯滤波核
       self.GFH = nn.Conv2d(1, 1, kernel_size=(1,filter_size), padding=(0,filter_size//2))
       self.GFV = nn.Conv2d(1, 1, kernel_size=(filter_size,1), padding=(filter_size//2,0))

# 设置 w 的值为 高斯平滑核, b 的值为 0.0
       init_parameter(self.GFH, generated_filters, np.array([0.0]))
       init_parameter(self.GFV, generated_filters.T, np.array([0.0]))

def forward(self, img):
       img_r = img[:,0:1]  # 取出RGB三个通道的数据
       img_g = img[:,1:2]
       img_b = img[:,2:3]

# 对图片的三个通道进行水平、垂直滤波
       blurred_img_r = self.GFV(self.GFH(img_r))
       blurred_img_g = self.GFV(self.GFH(img_g))
       blurred_img_b = self.GFV(self.GFH(img_b))

# 合并成一张图
       blurred_img = torch.stack([blurred_img_r, blurred_img_g, blurred_img_b], dim=1)
       blurred_img = torch.stack([torch.squeeze(blurred_img)])

return blurred_img

进行高斯平滑（模糊）之后的图片较原图更为模糊如下图右侧银币所示

完整代码见：gaussian_smoothing

（二）Sobel算子计算梯度

PAI = 3.1415926

class SobelFilterNet(nn.Module):
   def __init__(self) -> None:
       super(SobelFilterNet, self).__init__()
       sobel_filter = np.array([[-1, 0, 1],
                                [-2, 0, 2],
                                [-1, 0, 1]])
       self.SFH = nn.Conv2d(1, 1, kernel_size=sobel_filter.shape, padding=sobel_filter.shape[0]//2)
       self.SFV = nn.Conv2d(1, 1, kernel_size=sobel_filter.shape, padding=sobel_filter.shape[0]//2)

init_parameter(self.SFH, sobel_filter, np.array([0.0]))
       init_parameter(self.SFV, sobel_filter.T, np.array([0.0]))

def forward(self, img):
       img_r = img[:,0:1]
       img_g = img[:,1:2]
       img_b = img[:,2:3]

# # SFH(V): sobel filter of horizontal(vertical) 水平（竖直）方向的Sobel滤波
       grad_r_x = self.SFH(img_r)  # 通道 R 的 x 方向梯度
       grad_r_y = self.SFV(img_r)
       grad_g_x = self.SFH(img_g)
       grad_g_y = self.SFV(img_g)
       grad_b_x = self.SFH(img_b)
       grad_b_y = self.SFV(img_b)

# 计算强度（magnitude） 和 方向（orientation）
       magnitude_r = torch.sqrt(grad_r_x**2 + grad_r_y**2) # Gr^2 = Grx^2 + Gry^2
       magnitude_g = torch.sqrt(grad_g_x**2 + grad_g_y**2)
       magnitude_b = torch.sqrt(grad_b_x**2 + grad_b_y**2)

grad_magnitude = magnitude_r + magnitude_g + magnitude_b

grad_y = grad_r_y + grad_g_y + grad_b_y
       grad_x = grad_r_x + grad_g_x + grad_b_x

# tanθ = grad_y / grad_x 转化为角度 （方向角）
       grad_orientation = (torch.atan2(grad_y, grad_x) * (180.0 / PAI))
       grad_orientation =  torch.round(grad_orientation / 45.0) * 45.0  # 转化为 45 的倍数

return grad_magnitude, grad_orientation

将梯度强度当作图片进行输出，得到右下图最右侧图片，可知硬币的边缘区域梯度值较大（越大越亮）

完整代码见：sobel_filter

（三）非极大化抑制

非极大化抑制（NMS）的过程为：

• 将梯度强度矩阵grad_magnitude的每一点都作为中心像素点，与其同向或者反向的两个相邻点（共有8个）的梯度强度进行比较。

• 若中心点的梯度小于这两个方向上的梯度，则点中心的的梯度值设为0

进过上面的两个步骤，可以用一个像素的宽度替代了梯度屋脊效应，同时保留了屋脊的梯度强度（最大的梯度）。

class NonMaxSupression(nn.Module):
   def __init__(self) -> None:
       super(NonMaxSupression, self).__init__()

all_orient_magnitude = np.stack([filter_0, filter_45, filter_90, filter_135, filter_180, filter_225, filter_270, filter_315])

'''
       directional_filter功能见下面详细说明
       '''
       self.directional_filter = nn.Conv2d(1, 8, kernel_size=filter_0.shape, padding=filter_0.shape[-1] // 2)

init_parameter(self.directional_filter, all_filters[:, None, ...], np.zeros(shape=(all_filters.shape[0],)))

def forward(self, grad_magnitude, grad_orientation):

all_orient_magnitude = self.directional_filter(grad_magnitude)     # 当前点梯度分别与其其他8个方向邻域点做差（相当于二阶梯度）

'''
               \ 3|2 /
                \ | /
           4     \|/    1
       -----------|------------
           5     /|\    8
                / | \
               / 6|7 \
       注: 各个区域都是45度
       '''

positive_orient = (grad_orientation / 45) ％ 8             # 设置正方向的类型，一共有八种不同类型的方向
       negative_orient = ((grad_orientation / 45) + 4) ％ 8       # +4 = 4 * 45 = 180 即旋转180度(如 1 -(+4)-> 5)

height = positive_orient.size()[2]                        # 得到图片的宽高
       width = positive_orient.size()[3]
       pixel_count = height * width                                # 计算图片所有的像素点数
       pixel_offset = torch.FloatTensor([range(pixel_count)])

position = (positive_orient.view(-1).data * pixel_count + pixel_offset).squeeze() # 角度 * 像素数 + 像素所在位置

# 拿到图像中所有点与其正向邻域点的梯度的梯度（当前点梯度 - 正向邻域点梯度，根据其值与0的大小判断当前点是不是邻域内最大的）
       channel_select_filtered_positive = all_orient_magnitude.view(-1)[position.long()].view(1, height, width)

position = (negative_orient.view(-1).data * pixel_count + pixel_offset).squeeze()

# 拿到图像中所有点与其反向邻域点的梯度的梯度
       channel_select_filtered_negative = all_orient_magnitude.view(-1)[position.long()].view(1, height, width)

# 组合成两个通道
       channel_select_filtered = torch.stack([channel_select_filtered_positive, channel_select_filtered_negative])

is_max = channel_select_filtered.min(dim=0)[0] > 0.0 # 如果min{当前梯度-正向点梯度, 当前梯度-反向点梯度} > 0，则当前梯度最大
       is_max = torch.unsqueeze(is_max, dim=0)

thin_edges = grad_magnitude.clone()
       thin_edges[is_max==0] = 0.0

return thin_edges

directional_filter的用处是什么？

# 输入
tensor([[[[1., 1., 1.],  
         [1., 1., 1.],  
         [1., 1., 1.]]]])
# 输出
tensor([[[[0., 0., 1.],
         [0., 0., 1.],
         [0., 0., 1.]],

[[0., 0., 1.],
         [0., 0., 1.],
         [1., 1., 1.]],

[[0., 0., 0.],
         [0., 0., 0.],
         [1., 1., 1.]],

[[1., 0., 0.],
         [1., 0., 0.],
         [1., 1., 1.]],

[[1., 0., 0.],
         [1., 0., 0.],
         [1., 0., 0.]],

[[1., 1., 1.],
         [1., 0., 0.],
         [1., 0., 0.]],

[[1., 1., 1.],
         [0., 0., 0.],
         [0., 0., 0.]],

[[1., 1., 1.],
         [0., 0., 1.],
         [0., 0., 1.]]]], grad_fn=<ThnnConv2DBackward0>)

可知其获取输入的八个方向的梯度值（在当前项目的代码中，为获取当前点梯度与其它8个方向梯度之差）

根据梯度的强度和方向，将方向分成8个类别（即对于每一点有八个可能方向），如上代码中 "米" 型图所示。

下面给出计算当前点正向邻域的相邻点的梯度强度的过程（反向同理）

梯度方向grad_orientation: 0, 1,, 2, 3, 4, 5, 6, 7 (共有8哥方向)

各方向梯度强度all_orient_magnitude: [[..方向0的梯度..], [..方向1的梯度..], ..., [..方向7的梯度..]]

故对于方向为 i 的点，其在梯度强度中的位置为 all_orient_magnitude[i][x, y]，将all_orient_magnitude变化为一维向量后，对应的位置为position = current_orient × pixel_count + pixel_offset，我们就可以根据这个位置信息拿到当前点与其正向邻域点梯度强度之差（同理也可以拿到反向的）。

以下为辅助图示：

最后效果如下右侧图所示（左侧为未进行最大化抑制的图）

完整代码见：nonmax_supression

（四）滞后边缘跟踪

我们思考后发现，到目前为止仍有如下几个问题：

• 如果图像中有噪声，可能会出现边缘无关的点（伪边）

• 边缘点时阴时明

所以最后我们就需要进行滞后边缘跟踪了，其步骤如下：

• 设定两个阈值（一高一低），将梯度强度小于低阈值的像素点的梯度强度设为0，得到图像A

• 将梯度强度小于高阈值的像素点的梯度强度设为0，得到图像B

我们知道由于A的阈值较低，故边缘保留较完整，连续性较好，但是伪边可能也较多，B正好与A相反。

据此我们设想以B为基础，A为补充，通过递归追踪来补全B中边缺失的像素点。

to_bw = lambda image: (image > 0.0).astype(float)

class HysteresisThresholding(nn.Module):
   def __init__(self, low_threshold=1.0, high_threshold=3.0) -> None:
       super(HysteresisThresholding, self).__init__()
       self.low_threshold = low_threshold
       self.high_threshold = high_threshold

def thresholding(self, low_thresh: torch.Tensor, high_thresh: torch.Tensor):
       died = torch.zeros_like(low_thresh).squeeze()
       low_thresh = low_thresh.squeeze()
       final_image = high_thresh.squeeze().clone()

height = final_image.shape[0] - 1
       width = final_image.shape[1] - 1

def connected(x, y, gap = 1):
           right = x + gap
           bottom = y + gap
           left = x - gap
           top = y - gap

if left < 0 or top < 0 or right >= width or bottom >= height:
               return False

return final_image[top, left] > 0  or final_image[top, x] > 0 or final_image[top, right] > 0 \
               or final_image[y, left] > 0 or final_image[y, right] > 0 \
               or final_image[bottom, left] > 0 or final_image[bottom, x] > 0 or final_image[bottom, right] > 0

# 先高再宽
       def trace(x:int, y:int):
           right = x + 1
           bottom = y + 1
           left = x - 1
           top = y - 1
           if left < 0 or top < 0 or right >= width or bottom >= height or died[y, x] or final_image[y, x] > 0:
               return

pass_high = final_image[y, x] > 0.0
           pass_low = low_thresh[y, x] > 0.0

died[y, x] = True

if pass_high:
               died[y, x] = False
           elif pass_low and not pass_high:
               if connected(x, y) or connected(x, y, 2): # 如果其他方向有连接
                   final_image[y, x] = low_thresh[y, x]
                   died[y, x] = False

# 往回
           if final_image[y, x] > 0.0: # 当前点有连接
               if low_thresh[top, left] > 0: trace(left, top)
               if low_thresh[top, x] > 0: trace(x, top)    
               if low_thresh[top, right] > 0: trace(right, top)
               if low_thresh[y, left] > 0: trace(left, y)
               if low_thresh[bottom, left] > 0: trace(left, bottom)

# 往下
           trace(right, y)
           trace(x, bottom)
           trace(right, bottom)

for i in range(width):
           for j in range(height):
               trace(i, j)

final_image = final_image.unsqueeze(dim=0).unsqueeze(dim=0)

return final_image

def forward(self, thin_edges, grad_magnitude, grad_orientation):
       low_thresholded: torch.Tensor = thin_edges.clone()
       low_thresholded[thin_edges<self.low_threshold] = 0.0

high_threshold: torch.Tensor = thin_edges.clone()
       high_threshold[thin_edges<self.high_threshold] = 0.0

final_thresholded = self.thresholding(low_thresholded, high_threshold)

return low_thresholded, high_threshold, final_thresholded

如下图为依次为低阈值、高阈值的效果图

如下为滞后边缘跟踪后的效果图

可知其相对上方左侧图，一些伪边被消除了，相对右侧图，细节更加的丰富。

完整代码见：hysteresis_thresholding

来源：https://juejin.cn/post/7055129413616156686