从**DataWhale**中学习的。
1.1 简介
在图像处理领域中,特征点又被称为兴趣点或者角点,它通常具有旋转不变性和光照不变性和视角不变性等优点,是图像的重要特征之一,常被应用到目标匹配、目标跟踪、三维重建等应用中。点特征主要指图像中的明显点,如突出的角点、边缘端点、极值点等等,用于点特征提取的算子称为兴趣点提取(检测)算子,常用的有Harris角点检测、FAST特征检测、SIFT特征检测及SURF特征检测。
本次学习记录的是较为常用而且较为基础的Harris角点检测算法,它的思想以及数学理论能够很好地帮助我们了解兴趣点检测的相关原理。
1.2 内容介绍
1.2.1 基础知识
1.角点
上面这张图中:
- 中间那个矩形表示一个平坦区域,在各方向移动,窗口内像素值没有变化;
- 最下面那个矩形表示一个边缘特征(Edges),如果沿着垂直方向移动(梯度方向),像素值会发生改变;如果沿着边缘移动(平行于边缘) ,像素值不会发生变化;
- 左上角那个矩形表示一个角(Corners),不管你把它朝哪个方向移动,像素值都会发生很大变化。
上面的图比较形象的解释什么是角。然后对应下面的三种情况:
所以,最右边的图是一个角点。
1.3 角点类型
下图展示了不同角点的类型,可以发现:如果使用一个滑动窗口以角点为中心在图像上滑动,存在朝多个方向上的移动会引起该区域的像素值发生很大变化的现象。
1.4 图像梯度
“像素值发生很大变化”这一现象可以用图像梯度进行描述。在图像局部内,图像梯度越大表示该局部内像素值变化越大(灰度的变化率越大)。而图像的梯度在数学上可用微分或者导数来表示。对于数字图像来说,相当于是二维离散函数求梯度
先来理解下差分:
使用差分来近似导数:
G
x
(
x
,
y
)
=
H
(
x
+
1
,
y
)
−
H
(
x
−
1
,
y
)
G_x(x,y)=H(x+1,y)-H(x-1,y)
Gx(x,y)=H(x+1,y)−H(x−1,y)
G
y
(
x
,
y
)
=
H
(
x
,
y
+
1
)
−
H
(
x
,
y
−
1
)
G_y(x,y)=H(x,y+1)-H(x,y-1)
Gy(x,y)=H(x,y+1)−H(x,y−1)
1.4 Harris角点检测算法原理
1 算法思想
算法的核心是利用局部窗口在图像上进行移动,判断灰度是否发生较大的变化。如果窗口内的灰度值(在梯度图上)都有较大的变化,那么这个窗口所在区域就存在角点。
这样就可以将 Harris 角点检测算法分为以下三步:
- 当窗口(局部区域)同时向 x (水平)和 y(垂直) 两个方向移动时,计算窗口内部的像素值变化量 E ( x , y ) E(x,y) E(x,y) ;
- 对于每个窗口,都计算其对应的一个角点响应函数 R R R;
- 然后对该函数进行阈值处理,如果 R > t h r e s h o l d R > threshold R>threshold,表示该窗口对应一个角点特征。
2 第一步 — 建立数学模型
第一步是通过建立数学模型,确定哪些窗口会引起较大的灰度值变化。
让一个窗口的中心位于灰度图像的一个位置
(
x
,
y
)
(x,y)
(x,y),这个位置的像素灰度值为
I
(
x
,
y
)
I(x,y)
I(x,y) ,如果这个窗口分别向
x
x
x 和
y
y
y 方向移动一个小的位移
u
u
u和
v
v
v,到一个新的位置
(
x
+
u
,
y
+
v
)
(x+u,y+v)
(x+u,y+v) ,这个位置的像素灰度值就是
I
(
x
+
u
,
y
+
v
)
I(x+u,y+v)
I(x+u,y+v) 。
$就是窗口移动引起的灰度值的变化值。
设 w ( x , y ) w(x,y) w(x,y)为位置 ( x , y ) (x,y) (x,y)处的窗口函数,表示窗口内各像素的权重,最简单的就是把窗口内所有像素的权重都设为1,即一个均值滤波核。
当然,也可以把 w ( x , y ) w(x,y) w(x,y)设定为以窗口中心为原点的高斯分布,即一个高斯核。如果窗口中心点像素是角点,那么窗口移动前后,中心点的灰度值变化非常强烈,所以该点权重系数应该设大一点,表示该点对灰度变化的贡献较大;而离窗口中心(角点)较远的点,这些点的灰度变化比较小,于是将权重系数设小一点,表示该点对灰度变化的贡献较小。
则窗口在各个方向上移动
(
u
,
v
)
(u,v)
(u,v)所造成的像素灰度值的变化量公式如下:
若窗口内是一个角点,则
E
(
u
,
v
)
E(u,v)
E(u,v)的计算结果将会很大。为了提高计算效率,对上述公式进行简化,利用泰勒级数展开来得到这个公式的近似形式。
首先来看看二维的泰勒展开式:
所以
I
(
x
+
u
,
y
+
v
)
I(x+u,y+v)
I(x+u,y+v) 为:
I
(
x
+
u
,
y
+
v
)
=
I
(
x
,
y
)
+
u
I
x
+
v
I
y
I(x+u,y+v)=I(x,y)+uI_x+vI_y
I(x+u,y+v)=I(x,y)+uIx+vIy
其中
I
x
I_x
Ix和
I
y
I_y
Iy是
I
I
I的微分(偏导):
I
x
=
∂
I
(
x
,
y
)
∂
x
I_x=\frac{\partial I(x,y)}{\partial x}
Ix=∂x∂I(x,y),
I
y
=
∂
I
(
x
,
y
)
∂
y
I_y=\frac{\partial I(x,y)}{\partial y}
Iy=∂y∂I(x,y)
将
I
(
x
+
u
,
y
+
v
)
=
I
(
x
,
y
)
+
u
I
x
+
v
I
y
I(x+u,y+v)=I(x,y)+uI_x+vI_y
I(x+u,y+v)=I(x,y)+uIx+vIy代入
E
(
u
,
v
)
E(u,v)
E(u,v)可得:
提出 u 和 v ,得到最终的近似形式:
其中矩阵M为:
3 第二步—角点响应函数R
现在我们已经得到 E ( u , v ) E(u,v) E(u,v)的最终形式,别忘了我们的目的是要找到会引起较大的灰度值变化的那些窗口。
灰度值变化的大小则取决于矩阵M,M为梯度的协方差矩阵。在实际应用中为了能够应用更好的编程,所以定义了角点响应函数R,通过判定R大小来判断像素是否为角点。
计算每个窗口对应的得分(角点响应函数R定义):
其中
d
e
t
(
M
)
=
λ
1
λ
2
det(M)=\lambda_1\lambda_2
det(M)=λ1λ2是矩阵的行列式,
t
r
a
c
e
(
M
)
=
λ
1
+
λ
2
trace(M)=\lambda_1+\lambda_2
trace(M)=λ1+λ2 是矩阵的迹。
λ 1 λ1 λ1 和 λ 2 λ2 λ2 是矩阵 M M M的特征值, k k k是一个经验常数,在范围 (0.04, 0.06) 之间。
R R R的值取决于 M M M的特征值,对于角点 ∣ R ∣ |R| ∣R∣很大,平坦的区域 ∣ R ∣ |R| ∣R∣很小,边缘的 R R R为负值。
4 第三步——角点判定
根据 R 的值,将这个窗口所在的区域划分为平面、边缘或角点。这样得到了很多角点,为了得到最优的角点,我们还可以使用非极大值抑制(NMS)。 可以参考博客
注意:Harris 检测器具有旋转不变性,但不具有尺度不变性,也就是说尺度变化可能会导致角点变为边缘。想要尺度不变特性的话,可以关注SIFT特征。
因为特征值 λ1 和 λ2 决定了 R 的值,所以我们可以用特征值来决定一个窗口是平面、边缘还是角点:
- 平面::该窗口在平坦区域上滑动,窗口内的灰度值基本不会发生变化,所以 ∣ R ∣ |R| ∣R∣ 值非常小,在水平和竖直方向的变化量均较小,即 I x I_x Ix和 I y I_y Iy都较小,那么 λ1 和 λ2 都较小;
- 边缘: ∣ R ∣ |R| ∣R∣值为负数,仅在水平或竖直方向有较大的变化量,即 I x I_x Ix和 I y I_y Iy只有一个较大,也就是 λ1>>λ2 或 λ2>>λ1;
- 角点:[公式] 值很大,在水平、竖直两个方向上变化均较大的点,即 I x I_x Ix和 I y I_y Iy 都较大,也就是 λ1 和 λ2 都很大。
Harris 角点检测的结果是带有这些分数 R 的灰度图像,设定一个阈值,分数大于这个阈值的像素就对应角点。
5 补充一些
旋转不变性:Harris角点检测算子使用的是角点附近的区域灰度二阶矩矩阵。而二阶矩矩阵可以表示成一个椭圆,椭圆的长短轴正是二阶矩矩阵特征值平方根的倒数。当特征椭圆转动时,特征值并不发生变化,所以判断角点响应值R也不发生变化,由此说明Harris角点检测算子具有旋转不变性。
不具有尺度不变性: 如下图所示,当右图被缩小时,在检测窗口尺寸不变的前提下,在窗口内所包含图像的内容是完全不同的。左侧的图像可能被检测为边缘或曲线,而右侧的图像则可能被检测为一个角点。
1.5 基于opencv的实现
函数原型:cv2.cornerHarris(src, blockSize, ksize, k[, dst[, borderType]])
参数解释:
- src - 输入灰度图像,float32类型
- blockSize - 用于角点检测的邻域大小,就是上面提到的窗口的尺寸
- ksize - 用于计算梯度图的Sobel算子的尺寸
- k - 用于计算角点响应函数的参数k,取值范围常在0.04~0.06之间
代码示例:
import cv2 as cv
from matplotlib import pyplot as plt
import numpy as np
# detector parameters
block_size = 3
sobel_size = 3
k = 0.06
image = cv.imread('F:/DataWhale/5.jpg')
cop = image
print(image.shape)
# 获取图像的高度、宽度、通道数
height = image.shape[0]
width = image.shape[1]
channels = image.shape[2]
print("width: %s height: %s channels: %s"%(width, height, channels))
gray_img = cv.cvtColor(image, cv.COLOR_BGR2GRAY)
# modify the data type setting to 32-bit floating point
gray_img = np.float32(gray_img)
# detect the corners with appropriate values as input parameters
corners_img = cv.cornerHarris(gray_img, block_size, sobel_size, k)
# result is dilated for marking the corners, not necessary
kernel = cv.getStructuringElement(cv.MORPH_RECT,(3, 3))
dst = cv.dilate(corners_img, kernel)
# Threshold for an optimal value, marking the corners in Green
#image[corners_img>0.01*corners_img.max()] = [0,0,255]
for r in range(height):
for c in range(width):
pix=dst[r,c]
if pix>0.05*dst.max():
cv.circle(image,(c,r),2,(0,0,255),0)
image = cv.cvtColor(image, cv.COLOR_BGR2RGB)
plt.imshow(image)
plt.show()
结果:
参考的博客:
https://www.cnblogs.com/polly333/p/5416172.html
https://zhuanlan.zhihu.com/p/83064609
https://www.cnblogs.com/ronny/p/4009425.html