Harris角点检测是基于Moravec角点检测之上的, Moravec角点检测算子的思想其实特别简单,在图像上取一个W*W的“滑动窗口”,不断的移动这个窗口并检测窗口中的像素变化情况E。像素变化情况E可简单分为以下三种:A 如果在窗口中的图像是什么平坦的,那么E的变化不大。B 如果在窗口中的图像是一条边,那么在沿这条边滑动时E变化不大,而在沿垂直于这条边的方向滑动窗口时,E的变化会很大。 C 如果在窗口中的图像是一个角点时,窗口沿任何方向移动E的值都会发生很大变化。
用数学表达式就是:
(u,v)就表示四个移动方向(1,0)(1,1)(0,1)(-1,1),其中w为窗函数(window function),I为图像梯度,那么E就表示了灰度变化的剧烈程度
1977年,Moravec最先提出了如下的角点检测方法:
- 对于原始图像,取偏移量(Δx,Δy)为(1,0),(1,1),(0,1),(-1,1),分别计算每一像素点(xi,yi)的灰度变化
- 对于每一像素点(xi,yi),计算角点响应函数R(xi,yi)=min E
- 设定阈值T,将角点响应函数R(xi,yi)中低于T的值设为0
- 在窗口范围内进行非极大值抑制:遍历角点响应函数,若某个像素的角点响应函数在窗口内不是最大,该像素置0
- 选择非零点作为角点检测结果
Moravec角点检测的缺点
- 二值的窗口函数导致角点响应函数不够光滑
- 只在四个方向上计算灰度值变化,导致角点响应函数在多处都有较大响应
- 对于每个点只考虑E的最小值,导致算法对边缘有很强的反应
1988年,Harris和Plessey对Moravec的方法进行了改进,提出了经典的Harris角点检测算法。Harris首先将Moravec算法中的窗口函数由阶跃函数改为二维高斯函数,并通过泰勒展开考察微小移动,也就是说,如果要求E的最大值以明确角点,就可以令 u , v → 0 u,v \rightarrow 0 u,v→0,对E做泰勒展开,得
E
(
u
,
v
)
=
(
u
,
v
)
M
(
u
v
)
E(u, v)=(u, v) M \left( \begin{array}{l}{u} \\ {v}\end{array}\right)
E(u,v)=(u,v)M(uv)
M
=
∑
(
x
,
y
)
w
(
x
,
y
)
(
I
X
2
I
X
I
Y
I
X
I
Y
I
Y
2
)
=
(
∑
W
I
X
2
∑
W
I
X
I
Y
∑
W
I
X
I
Y
∑
W
I
Y
2
)
M=\sum_{(x, y)} w(x, y) \left( \begin{array}{cc}{I_{X}^{2}} & {I_{X} I_{Y}} \\ {I_{X} I_{Y}} & {I_{Y}^{2}}\end{array}\right)=\left( \begin{array}{cc}{\sum_{W} I_{X}^{2}} & {\sum_{W} I_{X} I_{Y}} \\ {\sum_{W} I_{X} I_{Y}} & {\sum_{W} I_{Y}^{2}}\end{array}\right)
M=(x,y)∑w(x,y)(IX2IXIYIXIYIY2)=(∑WIX2∑WIXIY∑WIXIY∑WIY2)
记
M
=
(
A
B
B
C
)
M=\left( \begin{array}{ll}{A} & {B} \\ {B} & {C}\end{array}\right)
M=(ABBC),则上式可以写成
A
u
2
+
2
B
u
v
+
C
v
2
=
E
A u^{2}+2 B u v+C v^{2}=E
Au2+2Buv+Cv2=E 的形式,这表示了一个椭圆,自相关矩阵M描述了图像局部区域的灰度变化趋势,可以通过椭圆的形状来判定角点。
下面解释下上面的内容:
1.窗口函数的两种形式:
2.泰勒公式展开:
,任何一个函数表达式,均可有泰勒公式进行展开,以逼近原函数,我们可以对下面函数进行一阶展开:
f
(
x
+
u
,
y
+
v
)
≈
f
(
x
,
y
)
+
u
f
x
(
x
,
y
)
+
v
f
y
(
x
,
y
)
f(x+u, y+v) \approx f(x, y)+u f_{x}(x, y)+v f_{y}(x, y)
f(x+u,y+v)≈f(x,y)+ufx(x,y)+vfy(x,y)
那么:
∑
[
I
(
x
+
u
,
y
+
v
)
−
I
(
x
,
y
)
]
2
\sum[I(x+u, y+v)-I(x, y)]^{2}
∑[I(x+u,y+v)−I(x,y)]2
≈
∑
[
I
(
x
,
y
)
+
u
I
x
+
v
I
y
−
I
(
x
,
y
)
]
2
\approx \sum\left[I(x, y)+u I_{x}+v I_{y}-I(x, y)\right]^{2}
≈∑[I(x,y)+uIx+vIy−I(x,y)]2
=
∑
u
2
I
x
2
+
2
u
v
I
x
I
y
+
v
2
I
y
2
=\sum u^{2} I_{x}^{2}+2 u v I_{x} I_{y}+v^{2} I_{y}^{2}
=∑u2Ix2+2uvIxIy+v2Iy2
=
∑
[
u
v
]
[
I
x
2
I
x
I
y
I
x
I
y
I
y
2
]
[
u
v
]
=\sum \left[ \begin{array}{cc}{u} & {v}\end{array}\right] \left[ \begin{array}{cc}{I_{x}^{2}} & {I_{x} I_{y}} \\ {I_{x} I_{y}} & {I_{y}^{2}}\end{array}\right] \left[ \begin{array}{l}{u} \\ {v}\end{array}\right]
=∑[uv][Ix2IxIyIxIyIy2][uv]
=
[
u
v
]
(
∑
[
I
x
2
I
x
I
y
I
x
I
y
I
y
2
]
)
[
u
v
]
=\left[ \begin{array}{ll}{u} & {v}\end{array}\right]\left(\sum \left[ \begin{array}{cc}{I_{x}^{2}} & {I_{x} I_{y}} \\ {I_{x} I_{y}} & {I_{y}^{2}}\end{array}\right]\right) \left[ \begin{array}{l}{u} \\ {v}\end{array}\right]
=[uv](∑[Ix2IxIyIxIyIy2])[uv]
把中间看成矩阵M,则:
E
(
u
,
v
)
≅
[
u
,
v
]
M
[
u
v
]
E(u, v) \cong[u, v] M \left[ \begin{array}{l}{u} \\ {v}\end{array}\right]
E(u,v)≅[u,v]M[uv]
M = [ ∑ I x 2 ∑ I x I y ∑ I x I y ∑ I y 2 ] = [ λ 1 0 0 λ 2 ] M=\left[ \begin{array}{cc}{\sum I_{x}^{2}} & {\sum I_{x} I_{y}} \\ {\sum I_{x} I_{y}} & {\sum I_{y}^{2}}\end{array}\right]=\left[ \begin{array}{cc}{\lambda_{1}} & {0} \\ {0} & {\lambda_{2}}\end{array}\right] M=[∑Ix2∑IxIy∑IxIy∑Iy2]=[λ100λ2]
现在回到椭圆函数那里,至于椭圆函数的得来,查到的是通过微分思想化简而来,具体的还没深究,
对于椭圆
A
x
2
+
2
B
x
y
+
C
y
2
=
1
A x^{2}+2 B x y+C y^{2}=1
Ax2+2Bxy+Cy2=1 ,设其半长轴和半短轴分别为a,b,那么
1
a
,
1
b
\frac{1}{\sqrt{a}}, \frac{1}{\sqrt{b}}
a1,b1 是矩阵M的特征值。(特征值的含义见添加链接描述)
将上面的点用一个椭圆轮廓包围
得到M的特征值有什么用呢?若用奇异值分解的观点看这个问题。由于Jacobian矩阵 J = ( I X I Y ) J=\left( \begin{array}{l}{I_{X}} \\ {I_{Y}}\end{array}\right) J=(IXIY),故M=JJ^T,这样M的特征值开根号后就是J的奇异值,因此M的特征值就可以体现IX和IY的相对大小。
现在我们需要定义角点响应函数,进一步进行区分,令(一般k=0.04 ~ 0.06)
R
=
det
(
M
)
−
k
tr
2
(
M
)
R=\operatorname{det}(M)-k \operatorname{tr}^{2}(M)
R=det(M)−ktr2(M)
M
=
[
∑
I
x
2
∑
I
x
I
y
∑
I
x
I
y
∑
I
y
2
]
=
[
λ
1
0
0
λ
2
]
M=\left[ \begin{array}{cc}{\sum I_{x}^{2}} & {\sum I_{x} I_{y}} \\ {\sum I_{x} I_{y}} & {\sum I_{y}^{2}}\end{array}\right]=\left[ \begin{array}{cc}{\lambda_{1}} & {0} \\ {0} & {\lambda_{2}}\end{array}\right]
M=[∑Ix2∑IxIy∑IxIy∑Iy2]=[λ100λ2]
其中:
det
(
M
)
=
λ
1
λ
2
=
∑
W
I
X
2
⋅
∑
W
I
Y
2
−
(
∑
W
I
X
I
Y
)
2
\operatorname{det}(M)=\lambda_{1} \lambda_{2}=\sum_{W} I_{X}^{2} \cdot \sum_{W} I_{Y}^{2}-\left(\sum_{W} I_{X} I_{Y}\right)^{2}
det(M)=λ1λ2=∑WIX2⋅∑WIY2−(∑WIXIY)2
tr
(
M
)
=
λ
1
+
λ
2
=
∑
W
I
X
2
+
∑
W
I
Y
2
\operatorname{tr}(M)=\lambda_{1}+\lambda_{2}=\sum_{W} I_{X}^{2}+\sum_{W}I_{Y}^{2}
tr(M)=λ1+λ2=∑WIX2+∑WIY2
一般增大k的值,将减小角点响应值R,降低角点检测的灵性,减少被检测角点的数量;减小k值,将增大角点响应值R,增加角点检测的灵敏性,增加被检测角点的数量。
参考博客:
https://www.jianshu.com/p/6ae09c3d226f
https://blog.csdn.net/newthinker_wei/article/details/45603583
https://zhuanlan.zhihu.com/p/42490675
https://zhuanlan.zhihu.com/p/36382429
https://blog.csdn.net/fengye2two/article/details/79119736