机器视觉——角点检测

什么是角点检测

在几何学里，我们会看到各种各样的三角形、多边形等，它们都有一个显著的特征：包含了角点信息。比如在三角形里，我们有三个角；在矩形里，我们有四个角。我们将找到这些图像特征的过程称为特征提取 (Feature Extraction)，我们之前所接触的Canny边缘检测也是特征提取的一种，其提取的是边缘特征，而我们这次所需要的是提取角点特征。现在，我们有一张图片需要去处理，要求：将图中书本的四个角点提取出来，以方便后续做透视变换等操作。这应该怎么样处理？

我们先对图像进行预处理操作

通过前面一系列知识的学习，我们可以将该RGB图片转换成灰度图像（本人比较懒，不想写代码，就用ps做处理了）：

随后，我们对该图像进行二值化操作：

可以看到该图像周围还是有很多的噪声，于是我们对图像进行形态学操作+中值滤波：

这样，我们的图像预处理工作就完成了，接下来是提取角点信息。我们通常有两种方法：Harris角点检测与Shi-Tomasi角点检测

Harris角点检测

我们怎么样去提取一张图像里的角点信息呢？我们可以使用一个滑动窗口对图像进行操作，如下图所示：

当窗口在图像上滑动时，我们可以根据窗口所在的区域分为如图的三种情况：

• 平坦区域：窗口内像素值几乎没有变化
• 边缘区域：沿平行于边缘的方向移动，像素值不会变化
• 角点区域：不管朝哪个方向移动，像素值都会发生变化‘

所以说，引起窗口内较大像素值变化的地方就是我们想要找寻的角点。那么我们该怎么样去找寻这些地方呢？

我们让一个窗口的中心位于灰度图像的一个位置((x,y))，这个位置的像素灰度值为(I(x,y))，如果这个窗口分别向x和y方向移动一个小的位移(Delta x)和(Delta y)，到一个新的位置((x+Delta x, y+Delta y))，这个位置的像素灰度值就是(I(x+Delta x, y+Delta y))，那么我们所说的灰度变化值就是([I(x+Delta x, y+Delta y)-I(x,y)])。我们将上式做个平方处理，因为灰度的变化有正有负；然后再给每个像素乘上一个权值(w(x,y))，表示图像在((x,y))处像素点的权值，而这个(w(x,y))也就是我们说的窗口函数。最简单的，我们可以把窗口内的权重都设置为1，窗口大小一般设置为3×3以上（根据实际情况而定）；有时我们也会将窗口内的像素权值设为一个高斯分布（二元正态分布）等（因为当窗口在角点位置移动时，中心点像素变化最大，所以将中心点的权值设为最大，周围逐渐下降，也就是一个高斯分布）。

这样，我们可以构造出一个函数，用来描述窗口在图像上移动时，窗口内像素灰度值的变化，该函数如下：

仅仅通过这个式子，我们还看不出来这个函数有什么性质，所以我们对该函数进行化简。

对于上面的过程，大家大概率看不太懂。别急，且听我一一道来

在 ((*)) 式中，我们将(Delta x和)(Delta y)看做极小量，对(I(x+Delta x,y+Delta y))使用了二元函数的全微分公式进行展开（其实就是对二元函数进行泰勒展开并省去高阶小量）（如果不是很明白全微分是什么东西，可以参考我的一篇博客：”https://www.cnblogs.com/Asaka-QianXiang/p/17172567.html”，自认为是讲的比较直观清楚的了）而对(I(x+Delta x,y+Delta y))进行如此操作的意图则是将变量(Delta x,Delta y)给分离出来，这样就能够方便我们对函数(E(Delta x, Delta y))进行分析，其实((**))式的意图也是如此。

在((**))式中，式子((Delta x^2 I_x^2 + Delta y^2 I_y^2 + 2Delta x Delta y I_x I_y))是以(Delta x,Delta y)为变量的二次型，于是我们写出该二次型的矩阵，也就是((**))式的东西。而我们这么处理的原因，就是对变量进行分离（用术语来说就是“解耦”）。我们的(Delta x, Delta y)是变量，而(I_x,I_y)都相当于常量，进行上述分解后，我们只需要研究矩阵：

为了方便书写，我们将该矩阵命名为(M)

在矩阵(M)里，(I_x,I_y)都是图像(I)在((x,y))处沿(x,y)轴的梯度（也就是(I)对(x,y)的偏导数，通俗点说就是斜率）。我们再仔细观察矩阵(M)，可以发现它是一个实对称矩阵。我们对实对称矩阵(M)进行对角化处理（其实也是“解耦”，再具体点说，是将矩阵(M)变换到正交坐标系里），也就是((***))式。对矩阵(M)对角化处理后，我们得到(M)的两个特征值(lambda_1)和(lambda_2)，这两个特征值的分别反映了图像在(x,y)方向的梯度信息。

是不是有点晕了，不要紧，我们能看懂最后的((***))式即可。当然，不要忘了我们最根本的目的：找到窗口像素变化很大的地方，也就是(Delta x,Delta y)稍微变一变，我的(E(Delta x,Delta y))就能产生剧烈变化。那么我们怎么将这个剧烈变化体现在我们的式子上呢？显然，就是我们的(Delta x, Delta y)的系数很大。在最终化简出来的((***))式中，(Delta x, Delta y)的系数是(R^{-1}
begin{bmatrix}
lambda_1 & 0
0 & lambda_2
end{bmatrix}
R)（令其为(A)），也就是除了(Delta x, Delta y)以外且与其内积的部分。学过线性代数中实对称矩阵对角化知识的同学都知道，我们的矩阵(R)里面包含的都是一堆正交的单位向量，也就是说，我们想知道系数(A)的大小，看矩阵(R)是没有用的，因为单位向量的模长是1，所以我们只能去分析矩阵(begin{bmatrix}
lambda_1 & 0
0 & lambda_2
end{bmatrix})。

我们刚刚也说了，(lambda_1)和(lambda_2)分别反映了图像在(x,y)方向的梯度信息。怎么去理解呢？我们可以考虑一个特殊情况：图像只包含垂直(x)轴或垂直(y)轴的边缘信息，也就是说，我们的(I_x)，(I_y)中有一个为0（不妨令(I_y)为0），因此矩阵M里的元素就只剩(sum_{(x,y)subset W}w(x,y) I_x^2)了，其实也就是(lambda_1)。即(lambda_1)特征值反映了x轴方向的梯度信息，后同理即可。

这样，我们的(E(Delta x, Delta y))函数就分析的差不多了，那么我们该怎么去综合这些信息，进而在图像上进行评估并找出角点呢？

我们还是要回到(lambda_1)和(lambda_2)。我们知道它们反映了图像梯度信息，那么当图像分别包含平坦、边缘、角点信息时，我们的(lambda_1)和(lambda_2)是怎样的呢？

首先，考虑平坦信息。此时，不管从哪个方向看，梯度几乎都是0，因为非常平坦，所以此时，(lambda_1,lambda_2)都很小

其次，考虑边缘信息。由于我们的((***))式其实已经将图像所有方向的梯度信息给变换到了(lambda_1,lambda_2)组成的矩阵上，也就是一组沿(x,y)坐标轴的正交向量（有点难理解），所以我们只要考虑水平边缘或垂直边缘即可。水平边缘时，沿x的梯度变化很大，即(lambda_2)很小（接近0），(lambda_1)很大；垂直边缘时，同理，(lambda_1)很小（接近0），(lambda_2)很大。

最后，考虑角点信息。此时不管沿哪个方向，梯度变化都很大，所以我们的(lambda_1,lambda_2)都很大。

有了上面的规律还不够，我们还要将其反映到一个函数上，也就是需要构造一个函数(R(lambda_1,lambda_2))，当(lambda_1,lambda_2)处于不同状态（角点和非角点）时，函数(R)会有显著不同的取值。（注意，这里的(R)函数与上文的(R)矩阵只是符号相同，意义完全不同）

我们有以下函数：

其中，R称为角点相应函数，det表示矩阵的行列式（也就是特征值之积），trace表示矩阵的迹（也就是对角线元素之和，即特征值之和），k为一个经验常数，在范围((0.04, 0.06))之间

用(lambda_1,lambda_2)表示函数R，即：

很容易验证：

• 当(lambda_1,lambda_2)都很小（接近0），也就是平坦信息时，R也很小

• 当(lambda_1,lambda_2)有一个很大，另一个很小时（边缘信息），R会有一个不是很小的数（设为(t)）

• 当(lambda_1,lambda_2)都很大且近似相等时（角点信息），R会有一个比(t)大的多的值

也就是说，我们只要设置一个阈值(threshold)，然后将滑动窗口遍历整个图像，找到(R(lambda_1, lambda_2)>threshold)的位置((x,y))，那么这个((x,y))就是服务器托管网我们需要找的角点位置。

梳理一下，Harris角点检测的整个流程就是：

• 输入灰度图像(Img)
• 使用(Sobel)算子等梯度算子计算图像上每一个点((x,y))处的梯度(I_x, I_y)
• 对图像里的每一个像素点，构造矩阵(M)，里面需要的参数有：窗口函数(w(x,y))，窗口里所有像素点的梯度信息(I_x,I_y)
• 求解角点响应函数(R)，筛选函数值大于threshold的点作为角点
• 统计所有角点信息并输出

当然，OpenCV里有API函数供我们直接调用，因此我们只需要调用现有函数即可。

函数原型：

cv2.cornerHarris(src, dst, blockSize, ksize, k, borderType=BORDER_DEFAULT)

参数：

src：		输入的灰度图像
dst：		Harris检测器输出的矩阵，大小等于输入图像，矩阵里每一个位置的值为角点响应函数R的值
blockSize：	角点检测中指定区域的大小（也就是W窗口的长或宽）
ksize：		Sobel求梯度操作中Sobel算子的大小
k：			角点响应函数中的经验常数k

代码例程（python）：

import cv2
import numpy as np

# 图像预处理
img = cv2.imread('Example4.png', flags=cv2.IMREAD_GRAYSCALE)
img_bgr = cv2.imread('Example.jpg')
img = cv2.resize(img, None, fx=0.5, fy=0.5)
img_bgr = cv2.resize(img_bgr, None, fx=0.5, fy=0.5)
kernel = np.ones((3,3), dtype=np.uint8)
img_open = cv2.morphologyEx(img, cv2.MORPH_OPEN, kernel, iterations=3)
img_close = cv2.morphologyEx(img_open, cv2.MORPH_CLOSE, kernel, iterations=10)
img_erode = cv2.morphologyEx(img_close, cv2.MORPH_ERODE, kernel, iterations=4)
img_blur = cv2.medianBlur(img_erode, 39)
cv2.imshow('blur', img_blur)
# Harris角点检测
harris_detector = cv2.cornerHarris(img_blur, 2, 3, 0.02)  # Harris角点检测器
dst = cv2.dilate(harris_detector, kernel, iterations=3)  # 为了使检测出来的角点更容易看见，所以膨胀
thresh = 0.011 * d服务器托管网st.max()  # 设置阈值
img_bgr[dst > thresh] = [0, 0, 255]  # 将角点标记为红色
cv2.imshow('dst', img_bgr)

cv2.waitKey(0)
cv2.destroyAllWindows()

效果：

可以发现，下方我们检测出来了很多角点，而且距离都很近，这时候我们就可以考虑使用非极大值抑制(Non-Maximum Suppression)算法。

算法内容可以自己上网查阅相关资料，这里简述下我的代码流程：

• 建立一个列表pt，存储检测到的角点坐标、Harris检测器响应函数大小信息
• 对列表pt按照响应大小进行降序排列
• 定义一个函数，能够返回图像上两个点距离的平方
• 设置阈值，如果两点距离平方大于这个阈值，就将响应较小的点丢弃，筛选出响应较大的点进入列表selected_pt

完整代码如下：

import cv2
import numpy as np

# 图像预处理
img = cv2.imread('Example4.png', flags=cv2.IMREAD_GRAYSCALE)
img_bgr = cv2.imread('Example.jpg')
img = cv2.resize(img, None, fx=0.5, fy=0.5)
img_bgr = cv2.resize(img_bgr, None, fx=0.5, fy=0.5)
kernel = np.ones((3,3), dtype=np.uint8)
img_open = cv2.morphologyEx(img, cv2.MORPH_OPEN, kernel, iterations=3)
img_close = cv2.morphologyEx(img_open, cv2.MORPH_CLOSE, kernel, iterations=10)
img_erode = cv2.morphologyEx(img_close, cv2.MORPH_ERODE, kernel, iterations=4)
img_blur = cv2.medianBlur(img_erode, 39)
cv2.imshow('blur', img_blur)
# Harris角点检测
harris_detector = cv2.cornerHarris(img_blur, 2, 3, 0.02)  # Harris角点检测器
thresh = 0.011 * harris_detector.max()  # 设置阈值
# NMS非极大值抑制部分
pt = []
for h in range(harris_detector.shape[0]):
    for w in range(harris_detector.shape[1]):
        if harris_detector[h, w] > thresh:
            pt.append([h, w, harris_detector[h,w]])
# 对pt中的元素按元素的第三个值降序排列
pt = sorted(pt, key=lambda x:x[2], reverse=True)


# 获取两点欧式距离的平方的函数
def getDist(pt1, pt2):
    return (pt1[0] - pt2[0]) * (pt1[0] - pt2[0]) + (pt1[1] - pt2[1]) * (pt1[1] - pt2[1])


selected_pt = []
t = 500  # 阈值，两点距离小于该值的话就取角点响应最强的点
selected_pt.append(pt[0])
for p in pt[1:]:
    if getDist(p, pt[0]) >= t:
        selected_pt.append(p)

print(selected_pt)
# 绘图，显示角点
for p in selected_pt:
    cv2.circle(img_bgr, (p[1], p[0]), 3, (0,0,255), 2)
cv2.imshow('dst', img_bgr)

cv2.waitKey(0)
cv2.destroyAllWindows()

输出的角点位置略有偏差，大家可以通过修改代码自行优化

Shi-Tomasi角点检测

相比较于Harris角点检测，Shi-Tomasi检测器修改了Harris检测器中的角点响应函数R。

Shi-Tomasi 发现角点的稳定性其实和矩阵 M 的较小特征值有关，于是直接用较小的那个特征值作为分数。这样就不用调整k值了。如果矩阵M的两个特征值中较小的那一个大于设定的阈值，那么这个点是角点；如果两个特征值都小于阈值，那较小的特征值也必定小于阈值，那这个点就是平坦区域的点；如果其中一个特征值大于阈值而另外一个特征值小于阈值，那么这个点就是边缘点。所以我们只需要判断矩阵M的两个特征值中较小的那一个特征值是否是大于阈值，如果大于阈值这个点就是强角点。
即角点响应函数修改为：

OpenCV中，我们使用cv2.goodFeaturesToTrack()函数进行Shi-Tomasi角点检测。

函数原型如下：

corners = cv2.goodFeaturesToTrack(img, 
                        maxCorners,
                        qualityLevel,
                        minDistance,
                        mask,
                        blockSize,
                        gradientSize,
                        useHarrisDetector=False, 
                        k=0.04
                       )

参数：

corners：		输出的角点向量
img：			输入图像
maxCorners：		能输出的最大角点个数（若maxCorners

应用实例：

import cv2
import numpy as np


# Shi-Tomasi角点检测部分
# 图像预处理
img = cv2.imread('Example4.png', flags=cv2.IMREAD_GRAYSCALE)
img_bgr = cv2.imread('Example.jpg')
img = cv2.resize(img, None, fx=0.5, fy=0.5)
img_bgr = cv2.resize(img_bgr, None, fx=0.5, fy=0.5)
kernel = np.ones((3,3), dtype=np.uint8)
img_open = cv2.morphologyEx(img, cv2.MORPH_OPEN, kernel, iterations=3)
img_close = cv2.morphologyEx(img_open, cv2.MORPH_CLOSE, kernel, iterations=10)
img_erode = cv2.morphologyEx(img_close, cv2.MORPH_ERODE, kernel, iterations=4)
img_blur = cv2.medianBlur(img_erode, 39)
cv2.imshow('blur', img_blur)
#  shi-tomasi角点检测
corners = cv2.goodFeaturesToTrack(img_blur, 4, 0.01, 100, blockSize=10)
print(corners)
#  显示角点
for i in corners:
    x,y = i.ravel()
    cv2.circle(img_bgr,(int(x),int(y)),5,(0,0,255),-1)
cv2.imshow('dst', img_bgr)

cv2.waitKey(0)
cv2.destroyAllWindows()

输出图像：

使用说明

一般来说，我们使用cv2.goodFeaturesToTrack函数来进行角点检测（即Shi-Tomasi检测器）。在OpenCV官方教程中，该函数说明如下：

（https://docs.opencv.org/4.1.0/dd/d1a/group__imgproc__feature.html#ga1d6bb77486c8f92d79c8793ad995d541）

我也仅仅说明一下需要注意的几个参数：

• 返回值corners是一个numpy数组。比如说，在上文的例子中，我们将corners 给打印出来：

想要处理这个corners数组，我们可以对里面的元素进行遍历，比如：

for i in corners:
    ...

但请注意，我们的i此时还是一个形如 [[811. 275.]] 的嵌套数组。想要将里面的数字提取出来，我们可以使用python内置的ravel方法：

for i in corners:
    x,y = i.ravel()

这样，我们就将检测器提取出的角点坐标转移到x,y变量中了，后续处理也方便不少。

• qualityLevel：调节该参数能够让我们控制角点数量：

我们的角点响应函数(R=min(lambda_1,lambda_2))，而且我们输入图像的每一个像素点都对应着一个响应函数R的值。假如说，我们的图像有一个像素点，其响应函数R的值最大（设为(m)），qualityLevel设为(0.01)，那么在整个图像中，只要一个像素点对应函数R的值大于(0.01*m)，那么这个点就会被认为是角点。

• blockSize：就是我们在角点检测原理一节中提及的滑动窗口W的长或宽。

调参Trick

一般来说，按以下要点调参即可：

• 想要多少角点，就将maxCorner设多大（建议一开始设大一些，然后逐渐缩小）
• qualityLevel一般0.01左右，qualityLevel的作用相当于阈值
• 检测出的角点密集就试试增大minDistance
• 需要提取的角点比较平滑（接近弧状），blockSize就给大些（如10），角点比较尖锐，blockSize就给小点（如3）

服务器托管，北京服务器托管，服务器租用 http://www.fwqtg.net
机房租用，北京机房租用，IDC机房托管， http://www.fwqtg.net

相关推荐: JavaWeb——第七章 前端工程化_2

第七章 前端工程化_2 六、Vue3视图渲染技术 6.1 模版语法 6.1.1 插值表达式和文本渲染 6.1.2 Attribute属性渲染 6.1.3 事件的绑定 6.2 响应式基础 6.2.1 响应式需求案例 6.2.2 响应式实现关键字ref 6.2.3…