专栏首页caoqi95的记录日志Python 实现 Canny 边缘检测算法

Python 实现 Canny 边缘检测算法

已同步:博客 | 公众号 | 语雀

Canny 边缘检测算法由计算机科学家 John F. Canny 于 1986 年提出的。其不仅提供了算法,还带来了一套边缘检测的理论,分阶段的解释如何实现边缘检测。Canny 检测算法包含下面几个阶段:

  • 灰度化
  • 高斯模糊
  • 计算图片梯度幅值
  • 非极大值抑制
  • 双阈值选取

灰度化

灰度化实际上是一种降维的操作,可以减少计算。如果算法不进行色彩相关的识别的话,不灰度化,也可以直接进行后面的阶段。

# 灰度化
def gray(self, img_path):
    """
    计算公式:
    Gray(i,j) = [R(i,j) + G(i,j) + B(i,j)] / 3
    or :
    Gray(i,j) = 0.299 * R(i,j) + 0.587 * G(i,j) + 0.114 * B(i,j)
    """
    # 读取图片
    img = plt.imread(img_path)
    # BGR 转换成 RGB 格式
    img_rgb = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
    # 灰度化
    img_gray = np.dot(img_rgb[...,:3], [0.299, 0.587, 0.114])
        
    return img_gray

高斯模糊

在实际的图片中,都会包含噪声。但有时候,图片中的噪声会导致图片中边缘信息的消失。对此的解决方案就是使用高斯平滑来减少噪声,即进行高斯模糊操作。该操作是一种滤波操作,与高斯分布有关,下面是一个二维的高斯函数,其中 (x, y) 为坐标,σ 为标准差:

进行高斯滤波之前,需要先得到一个高斯滤波器(kernel)。如何得到一个高斯滤波器?其实就是将高斯函数离散化,将滤波器中对应的横纵坐标索引代入高斯函数,即可得到对应的值。不同尺寸的滤波器,得到的值也不同,下面是 (2k+1)x(2k+1) 滤波器的计算公式 :

常用尺寸为 5x5,σ=1.4 的高斯滤波器。下面是 5x5 高斯滤波器的实现代码:

# 去除噪音 - 使用 5x5 的高斯滤波器
def smooth(self, img_gray):
        
    # 生成高斯滤波器
    """
    要生成一个 (2k+1)x(2k+1) 的高斯滤波器,滤波器的各个元素计算公式如下:
        
    H[i, j] = (1/(2*pi*sigma**2))*exp(-1/2*sigma**2((i-k-1)**2 + (j-k-1)**2))
    """
    sigma1 = sigma2 = 1.4
    gau_sum = 0
    gaussian = np.zeros([5, 5])
    for i in range(5):
        for j in range(5):
            gaussian[i, j] = math.exp((-1/(2*sigma1*sigma2))*(np.square(i-3) 
                                + np.square(j-3)))/(2*math.pi*sigma1*sigma2)
            gau_sum =  gau_sum + gaussian[i, j]
                
    # 归一化处理
    gaussian = gaussian / gau_sum
        
    # 高斯滤波
    W, H = img_gray.shape
    new_gray = np.zeros([W-5, H-5])
        
    for i in range(W-5):
        for j in range(H-5):
            new_gray[i, j] = np.sum(img_gray[i:i+5, j:j+5] * gaussian)
                
    return new_gray

图片梯度幅值

边缘是图像强度快速变化的地方,可以通过图像梯度幅值,即计算图像强度的一阶导数来识别这些地方。由于图片是离散的,可以用有限导数来近似图片的梯度:

图片梯度幅值为:

梯度方向为:

实现代码如下:

# 计算梯度幅值
def gradients(self, new_gray):
    """
    :type: image which after smooth
    :rtype: 
        dx: gradient in the x direction
        dy: gradient in the y direction
        M: gradient magnitude
        theta: gradient direction
    """
        
    W, H = new_gray.shape
    dx = np.zeros([W-1, H-1])
    dy = np.zeros([W-1, H-1])
    M = np.zeros([W-1, H-1])
    theta = np.zeros([W-1, H-1])
        
    for i in range(W-1):
        for j in range(H-1):
            dx[i, j] = new_gray[i+1, j] - new_gray[i, j]
            dy[i, j] = new_gray[i, j+1] - new_gray[i, j]
             # 图像梯度幅值作为图像强度值
            M[i, j] = np.sqrt(np.square(dx[i, j]) + np.square(dy[i, j]))
            # 计算  θ - artan(dx/dy)
            theta[i, j] = math.atan(dx[i, j] / (dy[i, j] + 0.000000001)) 
                
     return dx, dy, M, theta

非极大值抑制(NMS)

理想情况下,最终得到的边缘应该是很细的。因此,需要执行非极大值抑制以使边缘变细。原理很简单:遍历梯度矩阵上的所有点,并保留边缘方向上具有极大值的像素。

梯度方向与边缘方向相互垂直

下面说说 NMS 的细节内容。NMS 在 4 个方向上进行,分别是 0,90,45,135,没有角度包含两个领域,因此,一共用八个领域:上,下,左,右,左上,左下,右上,右下,如下图所示,C 周围的 8 个点就是其附近的八个领域。

这样做的好处是简单, 但是这种简化的方法无法达到最好的效果, 因为,自然图像中的边缘梯度方向不一定是沿着这四个方向的。因此,就有很大的必要进行插值,找出在一个像素点上最能吻合其所在梯度方向的两侧的像素值。

NMS 是要找出局部最大值,因此,需要将当前的像素的梯度,与其他方向进行比较。如下图所示,g1,g2,g3,g4 分别是 C 八个领域中的 4 个点,蓝线是 C 的梯度方向。如果 C 是局部最大值的话,C 点的梯度幅值就要大于梯度方向直线与 g1g2,g4g3 两个交点的梯度幅值,即大于点 dTemp1 和 dTemp2 的梯度幅值。上面提到这种方法无法达到最好的效果,因为 dTemp1 和 dTemp2 不是整像素,而是亚像素。亚像素的意思就是在两个物理像素之间还有像素。

那么,亚像素的梯度幅值怎么求?可以使用线性插值的方法,计算 dTemp1 在 g1,g2 之间的权重,就可以得到其梯度幅值。计算公式如下:

weight = |gx| / |gy| or |gy| / |gx|
dTemp1 = weight*g1 + (1-weight)*g2
dTemp2 = weight*g3 + (1-weight)*g4

下面两幅图是 y 方向梯度值比较大的情况,即梯度方向靠近 y 轴。所以,g2 和 g4 在 C 的上下位置,此时 weight = |gy| / |gx| 。左边的图是 x,y 方向梯度符号相同的情况,右边是 x,y 方向梯度符号相反的情况。

对于左边的图来说,以 C 点为当前位置 - d[i, j] ,那么 g2 在 C 的前一行,g4 在 C 的后一行,所以位置坐标是:<br />g2 = d[i-1, j];g4 = d[i+1, j]。根据左图的位置关系可以得到:g1 = d[i-1, j-1];g3 = d[i+1, j+1]。

同理,根据右图的位置关系可以得到:g1 = d[i-1, j+1];g3 = d[i+1, j-1]。

下面两幅图是 x 方向梯度值比较大的情况,即梯度方向靠近 x 轴。所以,g2 和 g4 在 C 的左右位置,此时 weight = |gy| / |gx| 。左边的图是 x,y 方向梯度符号相同的情况,右边是 x,y 方向梯度符号相反的情况。

由上面可知,可以得到如下信息:g2 = d[i, j-1];g4 = d[i, j+1]; 左图:g1 = d[i+1, j-1];g3 = d[i-1, j+1]; 右图:g1 = d[i-1, j-1];g3 = d[i+1, j+1]。

下面的这两幅图,可能会带来理解帮助:

然后,根据以上信息,代码实现如下:

def NMS(self, M, dx, dy):
        
    d = np.copy(M)
    W, H = M.shape
    NMS = np.copy(d)
    NMS[0, :] = NMS[W-1, :] = NMS[:, 0] = NMS[:, H-1] = 0
        
    for i in range(1, W-1):
        for j in range(1, H-1):
                
            # 如果当前梯度为0,该点就不是边缘点
            if M[i, j] == 0:
                NMS[i, j] = 0
                    
            else:
                gradX = dx[i, j] # 当前点 x 方向导数
                gradY = dy[i, j] # 当前点 y 方向导数
                gradTemp = d[i, j] # 当前梯度点
                    
                # 如果 y 方向梯度值比较大,说明导数方向趋向于 y 分量
                if np.abs(gradY) > np.abs(gradX):
                    weight = np.abs(gradX) / np.abs(gradY) # 权重
                    grad2 = d[i-1, j]
                    grad4 = d[i+1, j]
                        
                    # 如果 x, y 方向导数符号一致
                    # 像素点位置关系
                    # g1 g2
                    #    c
                    #    g4 g3
                    if gradX * gradY > 0:
                        grad1 = d[i-1, j-1]
                        grad3 = d[i+1, j+1]
                        
                    # 如果 x,y 方向导数符号相反
                    # 像素点位置关系
                    #    g2 g1
                    #    c
                    # g3 g4
                    else:
                        grad1 = d[i-1, j+1]
                        grad3 = d[i+1, j-1]
                    
                # 如果 x 方向梯度值比较大
                else:
                    weight = np.abs(gradY) / np.abs(gradX)
                    grad2 = d[i, j-1]
                    grad4 = d[i, j+1]
                        
                    # 如果 x, y 方向导数符号一致
                    # 像素点位置关系
                    #      g3
                    # g2 c g4
                    # g1
                    if gradX * gradY > 0:

                        grad1 = d[i+1, j-1]
                        grad3 = d[i-1, j+1]
                        
                    # 如果 x,y 方向导数符号相反
                    # 像素点位置关系
                    # g1
                    # g2 c g4
                    #      g3
                    else:
                        grad1 = d[i-1, j-1]
                        grad3 = d[i+1, j+1]
                        
                # 利用 grad1-grad4 对梯度进行插值
                gradTemp1 = weight * grad1 + (1 - weight) * grad2
                gradTemp2 = weight * grad3 + (1 - weight) * grad4
                    
                # 当前像素的梯度是局部的最大值,可能是边缘点
                if gradTemp >= gradTemp1 and gradTemp >= gradTemp2:
                    NMS[i, j] = gradTemp
                        
                else:
                    # 不可能是边缘点
                    NMS[i, j] = 0
                        
    return NMS

双阈值选取

这个阶段决定哪些边缘是真正的边缘,哪些边缘不是真正的边缘。为此,需要设置两个阈值,minVal 和 maxVal。梯度大于 maxVal 的任何边缘肯定是真边缘,而 minVal 以下的边缘肯定是非边缘,因此被丢弃。位于这两个阈值之间的边缘会基于其连通性而分类为边缘或非边缘,如果它们连接到“可靠边缘”像素,则它们被视为边缘的一部分。否则,也会被丢弃。

代码如下所示:

def double_threshold(self, NMS):
        
    W, H = NMS.shape
    DT = np.zeros([W, H])
        
    # 定义高低阈值
    TL = 0.1 * np.max(NMS)
    TH = 0.3 * np.max(NMS)
        
    for i in range(1, W-1):
        for j in range(1, H-1):
           # 双阈值选取
            if (NMS[i, j] < TL):
                DT[i, j] = 0
                    
            elif (NMS[i, j] > TH):
                DT[i, j] = 1
                    
           # 连接
            elif (NMS[i-1, j-1:j+1] < TH).any() or (NMS[i+1, j-1:j+1].any()
                    or (NMS[i, [j-1, j+1]] < TH).any()):
                DT[i, j] = 1
                    
        
    return DT 

边缘检测结果

经过以上 5 个过程,可以得到如下结果:

将其与 OpenCV,skimage 算法进行对比:

我个人感觉 OpenCV 的结果是最好的,其次是 Skimage 的结果。自己的算法结果有些地方还是蛮粗糙的。

完整的代码可以参见:caoqi95/CV_Learning/edge-detection

参考文献

[1]. 图像处理基础(4):高斯滤波器详解 [2]. Canny Edge Detection Step by Step in Python — Computer Vision [3]. canny 算子python实现 [4]. OpenCV-Canny Edge Detection

P.S: 代码主要参考了文章 [3],修改了一些错误并写成了类的形式。

PP.S: 文章有错误欢迎指出,互相学习 :)

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • POJ 1113 Wall(思维 计算几何 数学)

    给出平面上n个点的坐标。你需要建一个围墙,把所有的点围在里面,且围墙距所有点的距离不小于l。求围墙的最小长度。 n小于等于10^5

    attack
  • BZOJ1278: 向量vector(计算几何 随机化乱搞)

    首先我们可以按极角排序。然后对\(y\)轴上方/下方的加起来分别求模长取个最大值。。

    attack
  • 以后有面试官问你「跳跃表」,你就把这篇文章扔给他

    假如我们要用某种数据结构来维护一组有序的int型数据的集合,并且希望这个数据结构在插入、删除、查找等操作上能够尽可能着快速,那么,你会用什么样的数据结构呢?

    帅地
  • POJ 3304 Segments(直线与线段相交)

    我们考虑如果所有线段的投影有重合部分,那么我们可以在重合部分上做一条垂线经过所有点

    attack
  • 洛谷P4555 [国家集训队]最长双回文串(manacher 线段树)

    我的做法比较naive。。首先manacher预处理出以每个位置为中心的回文串的长度。然后枚举一个中间位置,现在要考虑的就是能覆盖到i - 1的回文串中 中心最...

    attack
  • 洛谷P1742 最小圆覆盖(计算几何)

    再枚举一个点\(j\),如果该点在圆内继续,否则用\(i, j\)构造出的圆替换出之前的圆。

    attack
  • BZOJ2564: 集合的面积(闵可夫斯基和 凸包)

    首先我们先分别求出给出的两个多边形的凸包。合并的时候直接拿个双指针扫一下,每次选最凸的点就行了。

    attack
  • UVa12298(生成函数简单应用、快速傅里叶变换)

    I have a set of super poker cards, consisting of an infinite number of cards. For...

    ACM算法日常
  • 每天一道剑指offer-构建乘积数组

    给定一个数组A[0,1,...,n-1],请构建一个数组B[0,1,...,n-1],其中B中的元素B[i]=A[0]A[1]...A[i-1]A[i+1].....

    乔戈里
  • leetcode 41 | 缺失的第一个正数 解题思路和方法

    思路:这个题刚開始是没有思路的,难就难在O(n)时间内常数量空间,所以此题较为考察思维敏捷性。其解题核心思想是将数组的第i位存正数i+1。最后再遍历一次就可以。

    ACM算法日常

扫码关注云+社区

领取腾讯云代金券