边缘识别

前言

视觉中有一个非常重要的模块叫边缘识别，顾名思义，这是用来帮助计算机如何识别图像边缘的技术。在这一节中，我们将介绍几种计算机视觉中常用的算法，并讲解原理和代码。

你也许需要的前置知识：线性代数+多元微积分

前置知识

颜色空间

是的，一张图片上的一个像素点可以由一个$1\times 3$的数组组成，其中每个分量代表了一种不可被继续分解的颜色——红、绿、蓝。也就是我们一般说的RGB颜色。而每种颜色有255个亮度，换句话说，255表示纯色（最大亮度），而0表示黑色（关闭显示这种颜色），比如，橙色可以由R = 200,G=100,B = 50混合而成。而每种颜色都是这三种不可分解颜色混合而成的，因此，利用线性代数的知识，我们知道这三种颜色生成一个空间，并且是以RGB为基的。

我们可以简单的验证颜色空间是一个线性空间，并且是以(1,0,0),(0,1,0)和(0,0,1)为基。(计算机中是以8bit去存储颜色的，所以值在0-255，显示中是以颜色的比例去调的，所以我们应该用0-1的值去做验证)

图像梯度

众所周知,数学上的梯度是表示一个函数在某个点的变化率（斜率）通过该点的变化率，我们可以得知这个点和上一个点或者是下一个点的关系，比如费马点，即对一元函数$f(x)$的导数$f’(x)$满足$f’(x) = 0$的点。那么这样的点要么是极大值要么是极小值。因为对该点的变化率前后都是要么大于它的，要么小于他的。

因此，若这些梯度很多的时候，我们可以看到这些箭头逐渐的组成函数的图形，可以看到从点0上去又下来。现在，我们回到图形上，一个图形由x和y两个方向组成，即水平和垂直两个方向，而图像梯度反应的是像素值(灰度或强度)在空间中的变化率，定义为：

对二维灰度图$I(x,y)$。它的梯度定义为图形对空间坐标偏导数组成的向量

$\nabla I = \left(I_x',I_y'\right)$

其中$I_x’$表示水平方向(x轴)上的变化率，$I_y’$表示垂直方向(y轴)上的变化率。接着，利用平行四边形定理，我们可以得到一个新的向量$\nabla I = I_x’ + I_y’$，新向量永远指向函数增长最大的地方。利用这个特性，在灰度图上我们就可以清楚的看到哪里是颜色变化最厉害的地方。那么这个最厉害的地方就有很大的可能是边界。

比如，我们定义函数$I(x,y)$是一个关于RGB内颜色运算的函数，那么他是下面这个图

它的图形梯度就是从黑色指向白色区域的向量，其上的函数是对颜色的值求变化率得到的

它刚好就是两个图形重叠的区域产生的边界。

现在我们开始真正的内容，进入到边缘检测的部分

Sobel算子

sobel算子是由 Irwin Sobel的名字命名的，它是下述的两个矩阵：

其中$G_x$来得到水平方向的导数，而$G_y$得到垂直方向的导数。

连续的梯度定义是取极限，而图形不同，因为图形的表示是离散的，所以实际上我们通过离散卷积的形式来实现类似于求导的操作。

Sobel算子的图像梯度计算步骤如下：

对输入图像进行卷积操作，分别计算图像在水平$G_x$和$G_y$方向的梯度。令$A = \{a_{ij}\}_{3\times 3}$是$3\times 3$的相邻像素矩阵，那么$A$上有9个相邻的像素点，那么它的梯度如下计算

$I_x' = A \times G_x $ 和 $I_y' = A \times G_y $

通过卷积核我们可以得知，其实$G_x$就是用后一个像素的颜色值去减前一个像素的颜色值，若颜色相差不大，这一减得到的值就趋近于0，否则这个值将增大。会得到一个不同于黑色非常明显的像素点。另一个也是一样的道理，基于此，卷积完得到的值就是边缘（高亮）+ 图片内点的像素（黑色）。

而计算完9个像素后，sobel一般使用步长为1的滑动窗口计算下一个矩阵，也就是将$a_{21}$的9个相邻像素点计算完之后就从$a_{22}$取9个相邻像素点计算。

接着，我们有两种方法去计算它的梯度的值，一种是欧式距离，即$G = \sqrt{I_x^2+I^2_y}$，另一种是出租车距离$G = \mid I_{x_2} - I_{x_1}\mid + \mid I_{y_2} - I_{y_1}\mid $，前者精度较高，后者较快。梯度的值表示边缘的强度，而梯度的方向表示边缘在哪。

只需要将得到的G标准化处理，那么我们就能得到最终的边缘图像。

代码编写

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
import numpy as np
import matplotlib.pyplot as plt

from PIL import Image
#G_x和G_y卷积核定义
G_x = np.array(
    [[-1,0,1],
     [-2,0,2],
     [-1,0,1]]
)

G_y = np.array(
    [
        [-1,-2,-1],
        [0,0,0],
        [1,2,1]
    ]
)
#图像处理函数
def sobel_test_fun(image_pth):
    img = Image.open(image_pth).convert('L')
    img_array = np.array(img,dtype=np.float32)
    height,width = img_array.shape 
    plt.imshow(img_array,cmap='gray')

    #初始化梯度数组

    grandient_x = np.zeros_like(img_array)#x
    grandient_y = np.zeros_like(img_array)#y
    g_magnitude = np.zeros_like(img_array) #x+y

    #卷积处理
    for i in range(1,height-1):
        for j in range(1,width -1):
            region = img_array[i-1:i+2, j-1:j+2]# 选取9个像素点
            
            grandient_x[i,j] = np.sum(region * G_x)
            grandient_y[i,j] = np.sum(region * G_y)

            g_magnitude[i,j] = np.sqrt(grandient_x[i,j]**2 + grandient_y[i,j]**2) #使用欧氏距离

    #归一化到0-255范围
    gradient_x = np.clip(grandient_x, 0, 255)
    gradient_y = np.clip(grandient_y, 0, 255)
    gradient_magnitude = np.clip(g_magnitude, 0, 255)

    return img_array, grandient_x,grandient_y,g_magnitude

def display(original,gx,gy,gm):
    plt.figure(figsize=(16, 12))
    
    plt.subplot(221)
    plt.imshow(original, cmap='gray')
    plt.title('orginal img')
    plt.axis('off')
    
    plt.subplot(222)
    plt.imshow(gx, cmap='gray')
    plt.title('img of x (Gx)')
    plt.axis('off')
    
    plt.subplot(223)
    plt.imshow(gy, cmap='gray')
    plt.title('img of y (Gy)')
    plt.axis('off')
    
    plt.subplot(224)
    plt.imshow(gm, cmap='gray')
    plt.title('edge of img')
    plt.axis('off')
    
    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    image_path = "your img"
    original,gx,gy,gm = sobel_test_fun(image_path)
    display(original,gx,gy,gm)

原图：

经过sobel算子得到的边缘图：

canny边缘检测算法

sobel算子是一种不错的检测方法，但不太够，容易受到噪声影响。比如有一笔小划痕在中间，用卷积核之后该划痕也会被划入边缘的集合内，这是不行的。现在我们介绍另一种算法，它的核心目标是低误测率、良好的定位性和最少的响应次数。它有5个主要步骤，即

  • 噪声抑制
  • 计算梯度
  • 非极大值抑制
  • 双阈值检测
  • 边缘连接

噪声抑制有：均值、中值、高斯、双边、引导滤波算法，这节我们讲高斯滤波

高斯滤波

高斯滤波基于正态分布函数，它是一个二维的函数

$G(x,y) = \frac{1}{2\pi\sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}}$

其中$\sigma^2$是方差，这样的运算也许有点复杂，如果它是下面这种形式就好了：

$G(x,y) = G(x)*G(y) = \cfrac{1}{2\pi\sigma^2}(e^{-x^2/2\sigma^2} * e^{-y^2/2\sigma^2})$

就好了，这意味着我们可以通过两次滤波来完成一次二维的高斯滤波，但这是可以的

proof:

$G(x) * G(y) = \begin{array}{ccc} \cfrac{1}{2\pi\sigma^2}e^{-(x^2/2\sigma^2 + y^2/2\sigma^2)} & = & G(x,y) \end{array}$

高斯模糊能起作用主要依赖于正态分布函数的标准差，在中心的周围1个标准差的距离是样本点最多的地方，2个标准差就开始减少，直到3个标准差，在大就没有意义了，因为3个标准差内包含99.7%的样本点。所以我们可以用这个特性加权。因此，我们可以如下的创建一个高斯模糊核，他用于和图片上的像素矩阵平均运算。离中心越远，则利用高斯曲线的特性快速下降权重，反之越靠近中心权重就越大，对结果影响也越大

给定标准差$\sigma = 1$，我们开始计算权重。给定一个$3\times 3$的矩阵$A = \{a_{ij}\}_{3\times 3}$，和$\sigma = 1$，矩阵的中心是$a_{22}$，以此为点$(0,0)$，则矩阵的第一个元素$a_{11}$相对于$a_{22}$是左上角，坐标为$(-1,1)$，右下角的是$(1,-1)$，右上角为$(1,1)$，左下角是$(-1,-1)$，我们只需要把坐标带入高斯函数，比如$G(a_{22}) = G(0,0) = 1/(2\pi)e^{-1} \approx 0.1592 $，那么高斯核就是

$G = \left(\begin{array}{ccc} 0.0545 & 0.0965 & 0.0545\\ 0.0965 & 0.1592 & 0.0965\\ 0.0545 & 0.0965 & 0.0545 \end{array}\right)$

接着需要归一化，若矩阵中的元素加起来不是1，则直接除以总数即可。

老套路，将高斯核$G$和图片的矩阵相乘即可得到平滑后的图片$G * A$，其中$A$是图片上提取的$3\times 3$像素。我们可以使用如下函数来调用高斯模糊

1
2
import cv2
blurred = cv2.GaussianBlur(img,(3,3),0) #第一个参数是图片，第二是高斯核大小，第三是sigma值，即标准差，设置为0则程序会自动计算合适的标准差值

接着和之前一样，使用sobel算子算出梯度强度和方向

$G = \sqrt{G_x^2+G_y^2}$ 和 $\theta = \arctan(G_y/G_x)$

最后一步，我们做范围的切割，$\theta \in [157.5,180] \cup [0,22.5]$的范围内，被定义为$0$。在$[22.5,67.5]$的范围内，$\theta$定义为$45$度，在$[67.5,112.5]$范围内被定义为$90$度。在$[112.5,157.5]$内被定义为$135$度