OpenCV

图像特征-harris角点检测

角点定义

角点检测又称为特征点检测，是图像处理和计算机视觉中用来获取局部特征点的一类方法，广泛应用于运动检测、图像匹配、视频跟踪等领域。角点可以简单地定义为轮廓之间的交点，严格地定义是在两个主方向上的特征点，即在两个方向上灰度变化剧烈。通常具有以下特征：

• 角点附近的像素点不论在梯度方向上还是梯度幅值上都存在着较大的变化
• 对于某一场景，当视角发生变化时，其任具备稳定性质的特征

算法原理

角点检测的基本思想就是用固定窗口在图像上沿各个方向进行滑动，比较滑动前后窗口中像素点的灰度变化，如果在任意方向上滑动窗口内都存在较大的灰度变化，则认为该窗口中存在角点；如果任何方向上都不变化，则是均匀区域；如果灰度只在一个方向上变化，则可能是图像边缘。

对于图像$I(x, y)$（主要指的是卷积核部分的区域）,当在点$(x, y)$处平移$(\Delta x,\Delta y)$后的相似性:

$c(x, y;\Delta x, \Delta y) = \sum_{(u, v) \in W(x, y)}W(u, v)(I(u, v)-I(u + \Delta x, v + \Delta y))^2\newline$

基于泰勒展开，对图像$I(x, y)$在平移$(\Delta x, \Delta y)$后进行一阶近似

$I(u+\Delta x, v + \Delta y)=I(u, v)+I_x(u, v)\Delta x+I_y(u, v)\Delta y + O(\Delta x^2, \Delta y^2) \approx I(u, v)+I_x(u, v)\Delta x+I_y(u, v)\Delta y\newline$

其中$I_x, I_y$是$I(x, y)$的偏导数

近似可得:

$c(x, y;\Delta x, \Delta y)\approx \sum_w(I_x(u, v)\Delta x+I_y(u, v)\Delta y)^2=[\Delta x, \Delta y]M(x, y)\begin{bmatrix}\Delta x\newline\Delta y\end{bmatrix}\newline$

其中$M$

$$ M(x, y)=\sum_w\begin{bmatrix}I_x(x, y)^2\ \ I_x(x, y)I_y(x, y)\newline I_x(x, y)I_y(x, y)\ \ I_y(x, y)^2\end{bmatrix}=\begin{bmatrix}\sum_w I_x(x, y)^2\ \ \sum_w I_x(x, y)I_y(x, y)\newline \sum_w I_x(x, y)I_y(x, y)\ \ \sum_w I_y(x, y)^2\end{bmatrix}=\begin{bmatrix}A\ \ C\newline C\ \ B\end{bmatrix}\newline $$

化简可得$c(x, y;\Delta x, \Delta y)\approx A\Delta x^2+2C\Delta x\Delta y+B\Delta y^2\newline$

其中$A=\sum_w I_x^2, B=\sum_wI_y^2, C=\sum_wI_xI_y\newline$

二次项函数本质上就是一个椭圆函数，椭圆方程为：$[\Delta x, \Delta y]M(x, y)\begin{bmatrix}\Delta x\newline\Delta y\end{bmatrix}=1\newline$

$W(x, y)$是以点$(x, y)$为中心的窗口，既可以是函数，也可以是高斯加权函数

• 平坦区域：在各方向移动，窗口内像素值均没有太大变化
• 边缘特征，如果沿着水平方向移动(梯度方向)，像素值会发生跳变；如果沿着边缘移动(平行于边缘) ，像素值不会发生变化
• 角：不管你把它朝哪个方向移动，像素值都会发生很大变化

算法流程

1. 滤波、平滑，避免出现阶跃函数。
2. 计算图像梯度。
3. 计算每个像素位置的Harris矩阵M。
4. 计算每个像素位置的角点响应函数R。
5. 设置门限R，寻找响应函数的局部最大值(非最大抑制)

代码

1

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

参数说明： src： 输入图像，须为float32型的单通道8位图像 blockSize： 邻域大小，即滑动窗口的尺寸 ksize： Sobel求导中使用的窗口大小，注意必须为奇数 k： Harris角点检测方程中的自由参数,推荐范围为[0.04, 0.06]，基本0.04

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

import cv2
import numpy as np

img = cv2.imread("img path")
print(f"img shape = {img.shape}") # (512, 512, 3)
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)  # 转化成灰度图
gray = np.array(gray, np.float32)
dst = cv2.cornerHarris(gray, 2, 3, 0.04) 
print(f"dst shape = {dst.shape}") # (512, 512)

img[dst > 0.01 * dst.max()] = [0, 0, 255]  # 筛选出来符合值
cv2.imshow('dst', img)
cv2.waitKey()
cv2.destroyAllWindows()

以上转载部分来自于CSDN幼儿园刘大壮的文章《图像特征-Harris角点检测》，侵删

Scale Invariant Feature Transform(SIFT)

图像尺度空间

在一定的范围内，无论物体是大还是小，人眼都可以分辨出来，然而计算机要有相同的能力却很难，所以要让机器能够对物体在不同尺度下有一个统一的认知，就需要考虑图像在不同尺度下都存在的特点

尺度空间的获取通常使用高斯模糊来实现

$I(x, y)$是一个图像 $$ \begin{align} &L(x, y, \sigma)=G(x, y, \sigma) * I(x, y)\newline &G(x, y, \sigma)=\frac {1}{2\pi \sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}}\newline \end{align} $$ $\sigma$控制模糊的程度,越大模糊效果越强

高斯差分金字塔（Difference of Gaussian)

• 首先先对图像进行差分

$DOG$定义公式

$$ D(x, y, \sigma)=[G(x, y, k\sigma)-G(x, y, \sigma)] * I(x, y)=L(x, y, k\sigma)-L(x, y, \sigma)\newline $$

DOG空间极值检测

为了寻找尺度空间的极值点，每个像素点要和其图像域（同一尺度空间）和尺度域（相邻的尺度空间）的所有相邻点进行比较，当其大于（或者小于）所有相邻点时，该点就是极值点。比如如果核的大小是$3\times 3$,中间的检测点就要和所在图像的$3\times 3$领域8个像素点，以及其相邻的上下两层的$3\times 3$领域18个像素点，共26个像素点进行比较

关键点的精确定位

这些候选关键点是$DOG$空间的局部极值点，而且这些极值点均为离散的点，精确定位极值点的一种方法是，对尺度空间$DOG$函数进行曲线拟合，计算其极值点，从而实现关键点的精确定位

可以用泰勒级数进行展开,离散的导数用差分代替微分计算 $$ f(x)\approx f(0)+f^`(0)x+\frac{f^{``}(0)}{2}x^2\newline $$

消除边界响应

$$ Hessian\ matrix\rightarrow H(x, y)=\begin{bmatrix}&D_{xx}(x, y)\ D_{xy}(x, y)\newline&D_{xy}(x, y)\ D_{yy}(x, y)\newline\end{bmatrix}\newline $$

令$\alpha=\lambda_{max}$为最大的特征值，$\beta=\lambda_{min}$为最小的特征值 $$ \begin{align} &T_r(H)=D_{xx}+D_{yy}=\alpha+\beta\newline &Det(H)=D_{xx}D_{yy}-(D_{xy})^2=\alpha\beta\newline &\frac {T_r(H)^2}{Det(H)}=\frac {(\alpha+\beta)^2}{\alpha\beta}=\frac {(\gamma + 1)^2}{\gamma}\newline \end{align} $$ $Lowe$在论文中给出了$\gamma=10$，也就是说对于主曲率比值大于10的特征点将会被删除

特征点的主方向

每个点$L(x, y)$的梯度的模$m(x, y)$以及方向$\theta(x, y)$ $$ \begin{align} &m(x, y)=\sqrt{[L(x+1, y)-L(x-1, y)]^2+[L(x, y+1)-L(x,, y-1)]^2}\newline &\theta(x, y)=\arctan(\frac {L(x, y+1)-L(x,, y-1)}{L(x+1, y)-L(x-1, y)})\newline \end{align} $$ 每个特征点可以得到三个信息$((x, y), \sigma, \theta)$,即位置，尺度和方向，具有多个方向的关键点可以被复制成多份，然后将方向值分别赋给复制后的特征点，一个特征点产生了多个坐标，尺度相等，但是方向不同的特征点

生成特征描述

在完成关键点的梯度计算后，使用直方图统计邻域内像素和梯度的方向

为了保证特征矢量的旋转不变性，要以特征点为中心，在附近邻域内将坐标旋转$\theta$角度，即将坐标轴旋转为特征点的主方向

$$ \begin{bmatrix} &x^’\newline &y^’\newline \end{bmatrix}=\begin{bmatrix} &\cos \theta\ -\sin \theta\newline &\sin \theta\ \cos \theta\newline \end{bmatrix}\begin{bmatrix} &x\newline &y\newline \end{bmatrix} $$

旋转之后的主方向为中心取$8\times 8$的窗口，求每个像素的梯度值和方向，箭头方向代表梯度方向，长度代表梯度幅值，然后利用高斯窗口对其进行加权运算，最后在每个$4\times 4$的小块上绘制8个方向的梯度直方图，计算每个梯度方向的累加值，即可形成一个种子点，即每个特征由四个种子点组成，每个种子点会有8个方向的向量信息

论文中建议对每个关键点使用$4\times 4$共16个种子点来描述，这样一个关键点就会产生128$[4(h),4(w),8(number\ of\ the\ vector)]$维的$SIFT$特征向量

OpenCV中的SIFT函数

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

import cv2

print(cv2.__version__)  # 4.6.0 not use the xfeatures2d.SIFT_create()
# ->(The patent has expired and the function has been deprecated), use the SIFT_create()

img = cv2.imread("img path")
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)  # 转化成灰度图
sift = cv2.SIFT_create()
kp = sift.detect(gray, None)  # 关键点

img = cv2.drawKeypoints(gray, kp, img)
cv2.imshow('drawKeyPoints', img)
cv2.waitKey()
cv2.destroyAllWindows()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

# 计算特征
kp, des = sift.compute(gray, kp)
print(np.array(kp).shape)
print(des.shape)
print(des[0])
'''
output
(789,)
(789, 128)
[ 47.  55.   9.  35. 117.  37.   0.   0. 143. 143.   2.   2.   2.   1.
   0.   0.  23.  23.   0.   0.   0.   0.   0.   0.   0.   0.   0.   0.
   0.   0.   0.   0.  36.  38.  64. 143.  40.   7.   1.   5. 143. 143.
  11.   8.   2.   0.   0.   9. 118.  23.   0.   0.   0.   0.   0.   1.
   1.   0.   0.   0.   0.   0.   0.   0.  36.   3.   9.  23.   7.  28.
  45.  68. 143.   8.   1.   0.   0.   0.   5. 143. 124.   1.   0.   0.
   0.   0.   0.  25.   1.   0.   0.   0.   0.   0.   0.   0.  52.   3.
   0.   0.   0.   4.  23. 119.  91.   0.   0.   0.   0.   0.   3. 143.
  26.   0.   0.   0.   0.   0.   1.  30.   0.   0.   0.   0.   0.   0.
   0.   0.]
'''

特征匹配

Brute-Force蛮力匹配

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

import cv2

img1 = cv2.imread("img1 path", 0)
img2 = cv2.imread("img2 path", 0)

sift = cv2.SIFT_create()

kp1, des1 = sift.detectAndCompute(img1, None)  # function -> 特征点，特征向量
kp2, des2 = sift.detectAndCompute(img2, None)

'''
crossCheck表示两个特征点要互相匹配，例如A中的第i个特征点与B中第j个特征点最近的
并且B中的第j个特征点到A中的第i个特征点也是(1对1)
此时计算方法为归一化数组(欧几里得距离) -> 默认参数:NORM_L2
'''

bf = cv2.BFMatcher_create(crossCheck=True)

# 一对一的匹配
matches = bf.match(des1, des2)
matches = sorted(matches, key=lambda x: x.distance)
img3 = cv2.drawMatches(img1, kp1, img2, kp2, matches[:50], None, flags=2)

K对最佳匹配

 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

import cv2


img1 = cv2.imread("img1 path", 0)
img2 = cv2.imread("img2 path", 0)

sift = cv2.SIFT_create()

kp1, des1 = sift.detectAndCompute(img1, None)  # function -> 特征点，特征向量
kp2, des2 = sift.detectAndCompute(img2, None)

bf = cv2.BFMatcher_create(crossCheck=False)

# k对匹配，注意上面的crossCheck要修正为False,否则1对1与下面k对匹配不符
matches = bf.knnMatch(des1, des2, k=2)
good = []

for m, n in matches:
    if m.distance < 0.75 * n.distance:
        good.append([m])

img3 = cv2.drawMatchesKnn(img1, kp1, img2, kp2, good, None, flags=2)
cv2.imshow("result", img3)
cv2.waitKey()
cv2.destroyAllWindows()
# 如果希望在数据结构方面优化，可以使用cv2.FlannBasedMatcher

缺点：会出现一些匹配错误的点

随机抽样一致算法($Random\ sample\ consensus,RANSAC$)

数据拟合

• 最小二乘法$\rightarrow$尽可能多的满足所有的数据点，部分异常数据（放在识别上就是异常点）会影响线的走向
• $RANSAC\rightarrow$能够减少异常点的干扰，得到的拟合线更符合实际

原理

• 选择初始样本点进行拟合，给定一个容忍范围，不断进行迭代
• 每一次拟合后，容差范围内都由对应的数据点数，找出数据点个数最多的情况，就是最终的拟合效果

单应性矩阵

$$ H=\begin{bmatrix} &h_{11}\ h_{12}\ h_{13}\newline &h_{21}\ h_{22}\ h_{23}\newline &h_{31}\ h_{32}\ 1\newline \end{bmatrix} $$

$$ \begin{bmatrix} &x^’\newline &y^’\newline &1 \end{bmatrix}=\begin{bmatrix} &h_{11}\ h_{12}\ h_{13}\newline &h_{21}\ h_{22}\ h_{23}\newline &h_{31}\ h_{32}\ 1\newline \end{bmatrix}\begin{bmatrix} &x\newline &y\newline &1\newline \end{bmatrix} $$

由线性代数知识可得，至少需要4对$(x^’,y^’)$的数据才能解该方程，为了防止误差点的干扰，筛选这些点需要使用$RANSAC$算法