Reference
https://blog.csdn.net/lbaihao/article/details/78919700
https://docs.opencv.org/4.x/dc/d0d/tutorial_py_features_harris.html
1 数学推导
Harris角点检测原理:
Harris算法利用窗口内图像灰度的自相关性进行角点判定,设置一个窗口,并在图像中移动,计算移动前与移动后窗口所在区域图像的自相关系数。
自相关函数计算如下,(x,y)为窗口中心位置,w(u,v)为权重(一般取高斯函数),L表示窗口,(u,v)表示窗口中的图像位置。
c
(
x
,
y
,
Δ
x
,
Δ
y
)
=
∑
(
u
,
v
)
∈
L
w
(
u
,
v
)
⋅
[
I
(
u
,
v
)
−
I
(
u
+
Δ
x
,
v
+
Δ
y
)
]
2
\mathrm{c}(\mathrm{x}, \mathrm{y}, \Delta \mathrm{x}, \Delta \mathrm{y})=\sum_{(\mathrm{u}, \mathrm{v}) \in \mathrm{L}} \mathrm{w}(\mathrm{u}, \mathrm{v}) \cdot[\mathrm{I}(\mathrm{u}, \mathrm{v})-\mathrm{I}(\mathrm{u}+\Delta \mathrm{x}, \mathrm{v}+\Delta \mathrm{y})]^{2}
c(x,y,Δx,Δy)=(u,v)∈L∑w(u,v)⋅[I(u,v)−I(u+Δx,v+Δy)]2
对其中的I二元函数进行泰勒展开:
I
(
u
+
Δ
x
,
v
+
Δ
y
)
=
I
(
u
,
v
)
+
I
x
(
u
,
v
)
⋅
Δ
x
+
I
y
(
u
,
v
)
⋅
Δ
y
+
⋯
+
R
n
(
u
,
v
)
≈
I
(
u
,
v
)
+
I
x
(
u
,
v
)
⋅
Δ
x
+
I
y
(
u
,
v
)
⋅
Δ
y
{\color{Brown} \begin{array}{c} \mathrm{I}(\mathrm{u}+\Delta \mathrm{x}, \mathrm{v}+\Delta \mathrm{y})=\mathrm{I}(\mathrm{u}, \mathrm{v})+\mathrm{I}_{\mathrm{x}}(\mathrm{u}, \mathrm{v}) \cdot \Delta \mathrm{x}+\mathrm{I}_{\mathrm{y}}(\mathrm{u}, \mathrm{v}) \cdot \Delta \mathrm{y}+\cdots+\mathrm{R}_{\mathrm{n}}(\mathrm{u}, \mathrm{v}) \\ \approx \mathrm{I}(\mathrm{u}, \mathrm{v})+\mathrm{I}_{\mathrm{x}}(\mathrm{u}, \mathrm{v}) \cdot \Delta \mathrm{x}+\mathrm{I}_{\mathrm{y}}(\mathrm{u}, \mathrm{v}) \cdot \Delta \mathrm{y} \end{array}}
I(u+Δx,v+Δy)=I(u,v)+Ix(u,v)⋅Δx+Iy(u,v)⋅Δy+⋯+Rn(u,v)≈I(u,v)+Ix(u,v)⋅Δx+Iy(u,v)⋅Δy
代入自相关函数,得:
c
(
x
,
y
,
Δ
x
,
Δ
y
)
=
∑
(
u
,
v
)
ϵ
L
w
(
u
,
v
)
⋅
[
I
(
u
,
v
)
−
I
(
u
+
Δ
x
,
v
+
Δ
y
)
]
2
≈
∑
(
u
,
v
)
∈
L
w
(
u
,
v
)
⋅
{
I
(
u
,
v
)
−
[
I
(
u
,
v
)
+
I
x
(
u
,
v
)
⋅
Δ
x
+
I
y
(
u
,
v
)
⋅
Δ
y
]
}
2
=
∑
(
u
,
v
)
∈
L
w
(
u
,
v
)
[
I
x
(
u
,
v
)
⋅
Δ
x
+
I
y
(
u
,
v
)
⋅
Δ
y
]
2
\begin{aligned} \mathrm{c}(\mathrm{x}, \mathrm{y}, \Delta \mathrm{x}, \Delta \mathrm{y}) &=\sum_{(\mathrm{u}, \mathrm{v}) \epsilon \mathrm{L}} \mathrm{w}(\mathrm{u}, \mathrm{v}) \cdot[\mathrm{I}(\mathrm{u}, \mathrm{v})-\mathrm{I}(\mathrm{u}+\Delta \mathrm{x}, \mathrm{v}+\Delta \mathrm{y})]^{2} \\ & \approx \sum_{(\mathrm{u}, \mathrm{v}) \in \mathrm{L}} \mathrm{w}(\mathrm{u}, \mathrm{v}) \cdot\left\{\mathrm{I}(\mathrm{u}, \mathrm{v})-\left[\mathrm{I}(\mathrm{u}, \mathrm{v})+\mathrm{I}_{\mathrm{x}}(\mathrm{u}, \mathrm{v}) \cdot \Delta \mathrm{x}+\mathrm{I}_{\mathrm{y}}(\mathrm{u}, \mathrm{v}) \cdot \Delta \mathrm{y}\right]\right\}^{2} \\ &=\sum_{(\mathrm{u}, \mathrm{v}) \in \mathrm{L}} \mathrm{w}(\mathrm{u}, \mathrm{v})\left[\mathrm{I}_{\mathrm{x}}(\mathrm{u}, \mathrm{v}) \cdot \Delta \mathrm{x}+\mathrm{I}_{\mathrm{y}}(\mathrm{u}, \mathrm{v}) \cdot \Delta \mathrm{y}\right]^{2} \end{aligned}
c(x,y,Δx,Δy)=(u,v)ϵL∑w(u,v)⋅[I(u,v)−I(u+Δx,v+Δy)]2≈(u,v)∈L∑w(u,v)⋅{I(u,v)−[I(u,v)+Ix(u,v)⋅Δx+Iy(u,v)⋅Δy]}2=(u,v)∈L∑w(u,v)[Ix(u,v)⋅Δx+Iy(u,v)⋅Δy]2
将平方项展开并写出矩阵形式,得:
[
I
x
(
u
,
v
)
⋅
Δ
x
+
I
y
(
u
,
v
)
⋅
Δ
y
]
2
=
I
2
x
(
u
,
v
)
Δ
x
2
+
I
2
y
(
u
,
v
)
Δ
y
2
+
2
I
x
(
u
,
v
)
I
y
(
u
,
v
)
Δ
x
Δ
y
=
I
2
x
(
u
,
v
)
Δ
x
2
+
I
2
y
(
u
,
v
)
Δ
y
2
+
2
I
x
(
u
,
v
)
I
y
(
u
,
v
)
Δ
x
Δ
y
=
[
I
x
2
(
u
,
v
)
Δ
x
+
I
x
(
u
,
v
)
I
y
(
u
,
v
)
Δ
y
I
y
2
(
u
,
v
)
Δ
y
+
I
x
(
u
,
v
)
I
y
(
u
,
v
)
Δ
x
]
⋅
[
Δ
x
Δ
y
]
=
[
Δ
x
Δ
y
]
⋅
[
I
x
2
(
u
,
v
)
I
x
(
u
,
v
)
I
y
(
u
,
v
)
I
x
(
u
,
v
)
I
y
(
u
,
v
)
I
y
2
(
u
,
v
)
]
⋅
[
Δ
x
Δ
y
]
\begin{array}{l} {\left[I_{x}(u, v) \cdot \Delta x+I_{y}(u, v) \cdot \Delta y\right]^{2}} \\ =I^{2}{ }_{x}(u, v) \Delta x^{2}+I^{2}{ }_{y}(u, v) \Delta y^{2}+2 I_{x}(u, v) I_{y}(u, v) \Delta x \Delta y \\ =I^{2}{ }_{x}(u, v) \Delta x^{2}+I^{2}{ }_{y}(u, v) \Delta y^{2}+2 I_{x}(u, v) I_{y}(u, v) \Delta x \Delta y \\ =\left[\begin{array}{lc} I_{x}^{2}(u, v) \Delta x+I_{x}(u, v) I_{y}(u, v) \Delta y & I_{y}^{2}(u, v) \Delta y+I_{x}(u, v) I_{y}(u, v) \Delta x \end{array}\right] \cdot\left[\begin{array}{l} \Delta x \\ \Delta y \end{array}\right] \\ =\left[\begin{array}{ll} \Delta x & \Delta y \end{array}\right] \cdot\left[\begin{array}{cc} I_{x}^{2}(u, v) & I_{x}(u, v) I_{y}(u, v) \\ I_{x}(u, v) I_{y}(u, v) & I_{y}^{2}(u, v) \end{array}\right] \cdot\left[\begin{array}{l} \Delta x \\ \Delta y \end{array}\right] \end{array}
[Ix(u,v)⋅Δx+Iy(u,v)⋅Δy]2=I2x(u,v)Δx2+I2y(u,v)Δy2+2Ix(u,v)Iy(u,v)ΔxΔy=I2x(u,v)Δx2+I2y(u,v)Δy2+2Ix(u,v)Iy(u,v)ΔxΔy=[Ix2(u,v)Δx+Ix(u,v)Iy(u,v)ΔyIy2(u,v)Δy+Ix(u,v)Iy(u,v)Δx]⋅[ΔxΔy]=[ΔxΔy]⋅[Ix2(u,v)Ix(u,v)Iy(u,v)Ix(u,v)Iy(u,v)Iy2(u,v)]⋅[ΔxΔy]
代入自相关函数,得:
c
(
x
,
y
,
Δ
x
,
Δ
y
)
≈
∑
(
u
,
v
)
∈
L
w
(
x
,
y
)
[
I
x
(
u
,
v
)
⋅
Δ
x
+
I
y
(
u
,
v
)
⋅
Δ
y
]
2
=
∑
(
u
,
v
)
∈
L
w
(
x
,
y
)
[
I
2
x
(
u
,
v
)
Δ
x
2
+
I
2
y
(
u
,
v
)
Δ
y
2
+
2
I
x
(
u
,
v
)
I
y
(
u
,
v
)
Δ
x
Δ
y
]
=
∑
(
u
,
v
)
∈
L
w
(
x
,
y
)
[
Δ
x
Δ
y
]
[
I
x
2
(
u
,
v
)
I
x
(
u
,
v
)
I
y
(
u
,
v
)
I
x
(
u
,
v
)
I
y
(
u
,
v
)
I
y
2
(
u
,
v
)
]
[
Δ
x
Δ
y
]
=
[
Δ
x
Δ
y
]
{
∑
(
u
,
v
)
∈
L
w
(
x
,
y
)
[
I
x
2
(
u
,
v
)
I
x
(
u
,
v
)
I
y
(
u
,
v
)
I
x
(
u
,
v
)
I
y
(
u
,
v
)
I
y
2
(
u
,
v
)
]
}
[
Δ
x
Δ
y
]
=
[
Δ
x
Δ
y
]
[
∑
(
u
,
v
)
ϵ
L
I
x
2
(
u
,
v
)
∑
(
u
,
v
)
ϵ
L
I
x
(
u
,
v
)
I
y
(
u
,
v
)
∑
(
u
,
v
)
∈
L
I
x
(
u
,
v
)
I
y
(
u
,
v
)
∑
(
u
,
v
)
ϵ
L
I
y
2
(
u
,
v
)
]
[
Δ
x
Δ
y
]
=
[
Δ
x
Δ
y
]
⋅
[
A
C
C
B
]
⋅
[
Δ
x
Δ
y
]
\begin{array}{l} c(x, y, \Delta x, \Delta y) \approx \sum_{(u, v) \in L} w(x, y)\left[I_{x}(u, v) \cdot \Delta x+I_{y}(u, v) \cdot \Delta y\right]^{2}\\ \\ =\sum_{(u, v) \in L} w(x, y)\left[I^{2}{ }_{x}(u, v) \Delta x^{2}+I^{2}{ }_{y}(u, v) \Delta y^{2}\right. \left.+2 I_{x}(u, v) I_{y}(u, v) \Delta x \Delta y\right]\\ \\ =\sum_{(u, v) \in L} w(x, y)\left[\begin{array}{ll} \Delta x & \Delta y \end{array}\right]\left[\begin{array}{cc} I_{x}^{2}(u, v) & I_{x}(u, v) I_{y}(u, v) \\ I_{x}(u, v) I_{y}(u, v) & I_{y}^{2}(u, v) \end{array}\right]\left[\begin{array}{l} \Delta x \\ \Delta y \end{array}\right]\\ \\ =\left[\begin{array}{ll} \Delta x & \Delta y \end{array}\right]\left\{\sum_{(u, v) \in L} w(x, y)\left[\begin{array}{cc} I_{x}^{2}(u, v) & I_{x}(u, v) I_{y}(u, v) \\ I_{x}(u, v) I_{y}(u, v) & I_{y}^{2}(u, v) \end{array}\right]\right\}\left[\begin{array}{l} \Delta x \\ \Delta y \end{array}\right]\\ \\ =\left[\begin{array}{ll} \Delta x & \Delta y \end{array}\right]\left[\begin{array}{cc} \sum_{(u, v) \epsilon L} I_{x}^{2}(u, v) & \sum_{(u, v) \epsilon L} I_{x}(u, v) I_{y}(u, v) \\ \sum_{(u, v) \in L} I_{x}(u, v) I_{y}(u, v) & \sum_{(u, v) \epsilon L} I_{y}^{2}(u, v) \end{array}\right]\left[\begin{array}{l} \Delta x \\ \Delta y \end{array}\right]\\ \\ =\left[\begin{array}{ll} \Delta x & \Delta y \end{array}\right] \cdot\left[\begin{array}{ll} A & C \\ C & B \end{array}\right] \cdot\left[\begin{array}{l} \Delta x \\ \Delta y \end{array}\right] \end{array}
c(x,y,Δx,Δy)≈∑(u,v)∈Lw(x,y)[Ix(u,v)⋅Δx+Iy(u,v)⋅Δy]2=∑(u,v)∈Lw(x,y)[I2x(u,v)Δx2+I2y(u,v)Δy2+2Ix(u,v)Iy(u,v)ΔxΔy]=∑(u,v)∈Lw(x,y)[ΔxΔy][Ix2(u,v)Ix(u,v)Iy(u,v)Ix(u,v)Iy(u,v)Iy2(u,v)][ΔxΔy]=[ΔxΔy]{∑(u,v)∈Lw(x,y)[Ix2(u,v)Ix(u,v)Iy(u,v)Ix(u,v)Iy(u,v)Iy2(u,v)]}[ΔxΔy]=[ΔxΔy][∑(u,v)ϵLIx2(u,v)∑(u,v)∈LIx(u,v)Iy(u,v)∑(u,v)ϵLIx(u,v)Iy(u,v)∑(u,v)ϵLIy2(u,v)][ΔxΔy]=[ΔxΔy]⋅[ACCB]⋅[ΔxΔy]
最终我们得到的公式为:
c
(
x
,
y
,
Δ
x
,
Δ
y
)
=
[
Δ
x
Δ
y
]
⋅
M
⋅
[
Δ
x
Δ
y
]
{\color{Red} \begin{array}{l} c(x, y, \Delta x, \Delta y) =\left[\begin{array}{ll} \Delta x & \Delta y \end{array}\right] \cdot M \cdot\left[\begin{array}{l} \Delta x \\ \Delta y \end{array}\right] \end{array}}
c(x,y,Δx,Δy)=[ΔxΔy]⋅M⋅[ΔxΔy]
M
=
∑
x
,
y
w
(
x
,
y
)
[
I
x
I
x
I
x
I
y
I
x
I
y
I
y
I
y
]
M = \sum_{x,y} w(x,y) \begin{bmatrix}I_x I_x & I_x I_y \\ I_x I_y & I_y I_y \end{bmatrix}
M=x,y∑w(x,y)[IxIxIxIyIxIyIyIy]
这里
I
x
I_x
Ix和
I
y
I_y
Iy是图像在x和y方向上的图像导数(使用cv2.Sobel()可以很容易得到它们)
椭圆的方程式
x
2
a
2
+
y
2
b
2
=
1
\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}=1
a2x2+b2y2=1
矩阵表达:
[
x
y
]
[
1
a
2
0
0
1
b
2
]
[
x
y
]
=
1
\left[\begin{array}{ll} x & y \end{array}\right]\left[\begin{array}{ll} \frac{1}{a^{2}} & 0 \\ 0 & \frac{1}{b^{2}} \end{array}\right]\left[\begin{array}{l} x \\ y \end{array}\right]=1
[xy][a2100b21][xy]=1
一个nxn的矩阵,可以求解其特征值,我们对上述中间矩阵求解特征值,可以得到
{
λ
1
=
1
a
2
λ
2
=
1
b
2
{
a
=
1
λ
1
b
=
1
λ
2
\left\{\begin{array} { l } { \lambda _ { 1 } = \frac { 1 } { a ^ { 2 } } } \\ { \lambda _ { 2 } = \frac { 1 } { b ^ { 2 } } } \end{array} \quad \left\{\begin{array}{l} a=\frac{1}{\sqrt{\lambda_{1}}} \\ b=\frac{1}{\sqrt{\lambda_{2}}} \end{array}\right.\right.
{λ1=a21λ2=b21{a=λ11b=λ21
特征值越大,轴半径越短。
M矩阵
c
(
x
,
y
,
Δ
x
,
Δ
y
)
=
[
Δ
x
Δ
y
]
⋅
M
⋅
[
Δ
x
Δ
y
]
{\color{Red} \begin{array}{l} c(x, y, \Delta x, \Delta y) =\left[\begin{array}{ll} \Delta x & \Delta y \end{array}\right] \cdot M \cdot\left[\begin{array}{l} \Delta x \\ \Delta y \end{array}\right] \end{array}}
c(x,y,Δx,Δy)=[ΔxΔy]⋅M⋅[ΔxΔy]
M
=
∑
x
,
y
w
(
x
,
y
)
[
I
x
I
x
I
x
I
y
I
x
I
y
I
y
I
y
]
M = \sum_{x,y} w(x,y) \begin{bmatrix}I_x I_x & I_x I_y \\ I_x I_y & I_y I_y \end{bmatrix}
M=x,y∑w(x,y)[IxIxIxIyIxIyIyIy]
我们发现自相关函数其实就是一个椭圆的矩阵表示。
那么,我们可以通过求解
M
M
M矩阵的特征值得到
λ
1
\lambda_1
λ1和
λ
2
\lambda_2
λ2。
- λ 1 \lambda_1 λ1和 λ 2 \lambda_2 λ2都特别小,轴长都很大,此时表示点在平坦区域。
- λ 1 \lambda_1 λ1和 λ 2 \lambda_2 λ2都特别大,而且比较接近,轴长都很小,可以认为是角点。
-
λ
1
\lambda_1
λ1和
λ
2
\lambda_2
λ2其中一个特别小,一个特别大,一个轴长,一个轴短,可以认为此时是边缘区域。
如果我们按照上面的方法,很明显要对每一个点都计算很多东西,速度太慢,因此提出了一个经验公式。
R
=
d
e
t
(
M
)
−
k
(
t
r
a
c
e
(
M
)
)
2
R=det(M)-k(trace(M))^2
R=det(M)−k(trace(M))2
其中,
- d e r ( M ) = λ 1 ⋅ λ 2 = A B − C 2 der(M)=\lambda_1\cdot\lambda_2=AB-C^2 der(M)=λ1⋅λ2=AB−C2
- t r a c e ( M ) = λ 1 + λ 2 = A + B trace(M)=\lambda_1+\lambda_2=A+B trace(M)=λ1+λ2=A+B
- λ 1 \lambda_1 λ1和 λ 2 \lambda_2 λ2是 M M M的特征值
- k k k是一个经验值,一般取0.02-0.04
2 Opencv 实现
2.1 函数 cv2.cornerHarris()【粗略了解就行】
2.1.1 参数
- img:输入图像,应该是float32的灰度图像
- blockSize: 角点检测所考虑的领域大小
- ksize: 使用Sobel导数的孔径参数
- k:Harris检测分数 R R R中 k k k的参数值
# -*-coding:utf-8-*-
import numpy as np
import cv2
filename = '1.jpg'
img = cv2.imread(filename)
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
# opencv需要输入是float32类型的
gray = np.float32(gray)
# 返回值是一个跟图片同样大小的矩阵分数图
dst = cv2.cornerHarris(gray, 5, 3, 0.04)
# 通过图像自身的分数图的大小关系,得到哪些像素点应该是角点,
# 这些像素点的值都变为[0,0,255],也就是红色。
img[dst>0.01*dst.max()]=[0,0,255]
cv2.imshow('dst', img)
cv2.waitKey(0)
2.2 cv2.goodFeaturesToTrack()
上面那个api的问题在于返回的不是角点的位置坐标,而是一个分数图,我们又直接用分数图处理了源图片。
cv2.goodFeaturesToTrack(),不仅支持Harris角点检测,也支持Shi Tomasi算法的角点检测。
2.2.1 参数
- gray_src:8位或者32位单通道灰度图像
- corners:位置点向量,保存的是检测到的角点的坐标
- max_corners:定义为可以检测到的角点的数量的最大值;
- qualityLevel:检测到的角点的质量等级,角点特征值小于qualityLevel最大特征值的点将被舍弃;
- minDistance:两角点间最小间距,以像素为单位;
- mask:指定检测区域,若检测整幅图像,则mask为空
- blockSize:计算自相关函数协方差矩阵的窗口大小;
- useHarrisDetector:是否使用Harris角点检测,为false,则使用Shi-Tomasi算子
- k:Harris角点检测的超参数,一般取0.02到0.04,userHarrisDetector为false,则不起作用。
# -*-coding:utf-8-*-
import numpy as np
import cv2
filename = '1.jpg'
img = cv2.imread(filename)
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
gray = np.float32(gray)
corners = cv2.goodFeaturesToTrack(gray, maxCorners=1000, qualityLevel=0.01, minDistance=10)
corners = np.int0(corners)
for i in corners:
x,y = i.ravel()
cv2.circle(img, (x, y), 2, (0, 0, 255), -1)
cv2.imshow('out.jpg', img)
cv2.waitKey(0)
亚像素
cv2.cornerSubPix()
参数:
- image:输入图像
- corners:输入角点的初始坐标
- winSize:搜索窗口边长的一半,如果winSize=(5,5),则一个大小为(52+1)(52+1)=1111的搜索窗口将被使用
- zeroZone:搜索区域中间的dead region边长的一半。如果值设置为(-1,-1)则表示没有这个区域
- cirteria:角点精准化迭代过程的终止条件,一般使用 c r i t e r i a = ( c v 2. T e r m C r i t e r i a E P S + c v 2. T e r m C r i t e r i a M A X I T E R , 100 , 0.001 ) criteria = (cv2.TermCriteria_EPS + cv2.TermCriteria_MAX_ITER, 100, 0.001) criteria=(cv2.TermCriteriaEPS+cv2.TermCriteriaMAXITER,100,0.001)