【Datawhale】计算机视觉下 —— Harris特征点检测

Harris特征点检测

openCV的安装

之前没有接触过openCV的小伙伴需要先在自己的环境下进行安装，因为笔者使用的是Mac系统和Anaconda环境，所以下面这个方案是面向Mac用户的。

# Mac系统中Anaconda下安装opencv
pip install opencv-python -i https://pypi.tuna.tsinghua.edu.cn/simple

# 检测是否安装成功
import cv2

基本概念

想要学习角点检测，首先应该知道什么是角点，它的应用场景，以及如何检测某个像素点是否为角点。

角点

严格的说，角点是物体边缘的拐点。在日常生活中，绝大多数的物体都是有棱角的，例如塔尖、三角形的顶角。那么这些有转折的，或者两个物体的交接处，都可以称作是角点。

但实际上，角点到目前为止都还没有明确的数学定义，所以我们只需要简单地认为角点是一个局部极值点，或者是二维图像亮度变化剧烈的点或图像边缘曲线上曲率的极大点。

应用

角点在如今已经广泛使用，如三维重建、图像匹配等工作。

如何判断角点

1. 以该像素点为中心，初始化一个滑动窗口（如3 ( imes) 3的滤波核），对它进行不同方向上的小范围移动（上下、左右、对角线），观察窗口内的像素灰度值变化情况。
2. 若无论怎么移动（滑动窗口内部始终包含此像素点），内部的像素值均无变化，或者变化很小，那么表明该点是一个平坦点（flat）。
3. 若此区域内的像素值在移动的过程中，在某个方向上的变化较大（单一方向），表明该点是个边缘点（edge）。
4. 若该区域内的像素值在移动的过程中，对多个方向上的像素值变化都很明显，则说明这个点是个角点（corner）。
5. 上述表述中的“明显”显然是我们主观的定义，那么到底多大的变化才算明显呢？我们需要对变化结果值经过一定的处理后得到最终结果(T)，并对它做一个阈值处理，假设此时的阈值为(K)，那么当(T > K)时，则该像素点为角点。

模型搭建与理解

前面我们已经从概念上理解了角点是什么，同时也明白了角点与普通像素点在邻域窗口中变化结果的不同。那么，到底应该具体化地表示移动下的窗口内部像素变化呢？这就涉及了数学建模的过程。

下面，我们给出在此邻域窗口内的像素灰度变化计算表示：

[E(u, v) = sum_{x, y} w(x, y) [I(x + u, y + v) - I(x, y)] ^ 2 ]

模型元素

已知(x, y)分别为该像素点的横纵坐标，(u, v)是窗口在水平、数值方向上的偏移。此时上述公式可拆解为四部分：

(E(u, v))：表示平移后的窗口内部像素变化值。假设此时的窗口是(3 imes 3)的滤波核，那么(E(u, v))就等于9个像素点的灰度变化之和。

(w(x, y))：窗口内各像素的权重。比较常用的是高斯核，此时考虑的是若该像素点为角点，那么它的移动一定使得灰度差值很大，所以该点的差值理应比其他点所占的权重更大。

(I(x, y))：函数(I)可求得某坐标处的像素点的灰度值。所以此时求得的是坐标为((x, y))时的灰度值。

(I(x+u, y+v) - I(x, y))：移动后对应像素的变化值。比如，此时该邻域窗口朝右下方移动，有(u = 1, v = 1)，那么对于中心点的灰度变化为 (I(x+1, y+1) - I(x, y))。

⚠️ 之所以在差值外加上平方，是因为在计算过程中，我们并不能保证所有的差值都是正数，可能会出现负值。但实际上我们并不关心它的正负性，所以加平方后的求和可以消除这一影响。

公式推导

首先要明确一个目的，就是我们想要寻找一个像素点，使得我们的(u, v)无论怎样取值，(E(u, v))都是变化比较大的，这个像素点就是我们要找的角点。显然现在的公式无法让我们直观地看出在什么条件下，可以使得该值尽可能地大，因此我们将使用泰勒公式对其进行一定的化简，然后变形。

将(I(x+u,y+v)=I(x,y)+uI_x+vI_y)带入上述公式，得到

[E(u, v) = sum _{x, y} w(x, y) imes (u^2I_x + v^2I_y + 2uvI_xI_y) ]

其中(I_x=frac{partial I(x,y)}{partial x})，(I_y=frac{partial I(x,y)}{partial y}) 。

提出公式内部的(u, v)后得到：

[E(u, v) = [u, v] sum _{x, y} w(x, y) egin{bmatrix} I_x^2 & I_xI_y \ I_xI_y & I_y^2 \ end{bmatrix} egin{pmatrix} u \ v \ end{pmatrix} ]

将求和公式移入矩阵内部，得到

[egin{bmatrix} sum _{x, y} w(x, y)I_x^2 & sum _{x, y} w(x, y)I_xI_y \ sum _{x, y} w(x, y)I_xI_y & sum _{x, y} w(x, y)I_y^2 \ end{bmatrix} = egin{bmatrix} A & C \ C & B \ end{bmatrix} = M ]

所以最后得到：

[E(u, v) = [u, v] M egin{pmatrix} u \ v \ end{pmatrix} ]

那么这时，成功将该模型的数值变化与(M)的变化相关联。在(u,v)值不变的情况下，(E(u, v))受(M)的影响较大。

在线性代数中，有概念：实对称矩阵可以正交相似对角化。那么根据公式4，我们知道(M)是一个实对称矩阵，所以有以下变形：

[M = Pegin{bmatrix} lambda_1 & 0\ 0 & lambda_2 \ end{bmatrix} P^{-1} ]

且其中(P)是正交矩阵，有性质如下:

[P^T = P^{-1} ]

综合公式5、6、7，可以得到

[E(u, v) = [u, v] Pegin{bmatrix} lambda_1 & 0\ 0 & lambda_2 \ end{bmatrix} P^T [u, v]^T = [u', v'] egin{bmatrix} lambda_1 & 0\ 0 & lambda_2 \ end{bmatrix}[u', v'] ^ T ]

[E(u, v) = lambda_1(u')^2 + lambda_2(v')^2 = frac{(u')^2}{frac{1}{lambda_1}} + frac{(v')^2}{frac{1}{lambda_2}} ]

所以，我们终于可以得到结果，(E(u, v))实际上就是一个椭圆呀，且根据咱们高中的知识，知道横轴是(lambda{min}^{-frac{1}{2}})，纵轴为(lambda{max}^{-frac{1}{2}})。如下图所示：

那么，前面提到我们的目标是找到一个像素点使得(E(u, v))最大，经过化简后，我们的目标转换成了找到一组(lambda_1, lambda_2)使得椭圆最大。

角点检测函数

这节就来讨论一下如何使得椭圆最大吧。从高中知识来看，显然当该椭圆的横纵轴同时很大时，它才会是最大的。那么，怎样的函数才能确保当且仅当(lambda_1, lambda_2)同时很大的情况下，函数值是最大的呢？

在这里我们定义下列公式：

[det(M) = lambda_1 lambda_2 ]

[trace (M) = lambda_1 + lambda_2 ]

其中det表示矩阵行列式，trace表示矩阵的迹。

随后，即可定义我们的角点检测函数为：

[T = det(M) - k(trace(M))^2 ]

在这个公式中，k值的设定是openCV中的预设，大概也经过大量实验之后得到的“经验”值。我们将借助这个角点检测函数，对一个像素点是否为角点做一个最终判定。

角点的最终判定

在实际的角点判定中，应该假设一个适当的阈值(K)。当我们求得的(T > K)时，该像素点将被认为是一个角点。

由于结构张量(M)是由(I_x，I_y)构成，它与特征值成正比关系，那么其特征值正好可以反映两者的情况。特征值 (λ_1) 和(λ_2)决定了 T 的值，所以我们可以用特征值来决定一个窗口是平面、边缘还是角点，如下图所示：

前面提到，实际上(E(u, v))的分布就像是一个椭圆，可能从这张图还不能清楚的看出，那么就借助下面这张实验结果图具体展示一下吧！

1. 右上角的图中(I_x，I_y)都比较小，也就是(λ_1) 和(λ_2)都较小，红色实线框住了所有的蓝点，发现它大体是一个面积较小的圆。

2. 右上角的图中(I_x，I_y)都比较大，也就是(λ_1) 和(λ_2)都较大，发现它是一个面积较大的、接近圆形的椭圆。

3. 右上角的图中(I_x > I_y)，也就是(λ_1 > λ_2)，发现它是一个扁平的椭圆，且横轴在X轴上。

那么现在，根据上述实验结果与分块图，可以总结如下结论，假设现在有阈值(K)：

1. 若(lambda_1、lambda_2)都很小，使得(0< T < K)，表示那么图像窗口在所有方向上的移动都无明显的灰度变化，所以该点被认为是处在平面上。

2. 若(lambda_1 >> lambda_2 || lambda_2 >> lambda_1)，使得(T < 1)，表明图像窗口仅在水平或竖直方向有较大的变化量，所以该点被认为是处在边缘处。

3. 若(lambda_1、lambda_2)都很大，使得(T > K)，表示那么图像窗口在各方向上的移动灰度变化明显，所以该点被认为是角点。

代码实现

导入环境包

在进行角点检测的过程中，我们需要使用的是以下几个环境包，主要负责图像的读取和显示。

import cv2
import matplotlib.pyplot as plt
import numpy as np

读取图片

我将希望进行检测的图片放在img包下，并使用cv2的cv2.imread()方法进行读取。

img = cv2.imread("../img/cv_1.jpg")
print(img.shape)  # (450, 720, 3)

灰度转换与角点检测函数

在opencv中有提供实现 Harris 角点检测的函数 cv2.cornerHarris，我们直接调用的就可以，非常方便。

函数原型：cv2.cornerHarris(src, blockSize, ksize, k[, dst[, borderType]])

对于每一个像素 (x,y)，在 (blockSize x blockSize) 邻域内，计算梯度图的协方差矩阵 (M(x,y))，然后通过上面第二步中的角点响应函数得到结果图。图像中的角点可以为该结果图的局部最大值。

即可以得到输出图中的局部最大值，这些值就对应图像中的角点。

参数解释：

• src - 输入灰度图像，float32类型
• blockSize - 用于角点检测的邻域大小，就是上面提到的窗口的尺寸
• ksize - 用于计算梯度图的Sobel算子的尺寸
• k - 用于计算角点响应函数的参数k，取值范围常在0.04~0.06之间

block_size = 2
sobel_size = 3
k = 0.04

gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) # 转化为灰度图
gray = np.float32(gray)
dst = cv2.cornerHarris(gray, block_size, sobel_size, k)

阈值设定

此次我们将阈值设置为图片中所有像素点中最大灰度值的0.01倍。假如某点灰度大于该值，则直接将颜色修改为[0, 0, 255]，为红色点。

img[dst > 0.01 * dst.max()] = [0, 0, 255]

结果显示

img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB) # 重新转化为彩色图
plt.imshow(img)
plt.show()

最后结果如下图所示（可能需要放大才能看清...）：

参考文献

https://blog.csdn.net/qq_40855100/article/details/106603051

https://blog.csdn.net/majinlei121/article/details/51043359