Harris角点检测原理
Harris角点检测是一种计算机视觉中常用的角点检测算法，用于在图像中检测出角点特征。角点通常被定义为两条边的交点，或者说，角点的局部邻域应该具有两个不同区域的不同方向的边界。Harris角点检测算法是最常用且最基础的角点检测算法之一。
Harris角点检测算法的原理是通过计算图像中每个像素点的响应函数来判断其是否为角点。该响应函数基于图像的灰度变化和局部窗口的协方差矩阵，通过计算特征值来确定角点的位置。当特征值较大时，表示该像素点处存在角点。
优点： 对图像旋转、尺度变化和亮度变化具有一定的不变性，且计算简单快速。
当一幅图像
f
f
f，检测窗口在图像上滑动，当检测窗口在点
(
x
,
y
)
(x,y)
(x,y)处位移
(
Δ
x
,
Δ
y
)
(\Delta x,\Delta y)
(Δx,Δy)后，它的自相关函数为：
C
(
x
,
y
)
=
∑
x
i
∈
w
∑
y
i
∈
w
(
f
(
x
i
,
y
i
)
−
f
(
x
i
−
Δ
x
,
y
i
−
Δ
y
)
)
2
C(x,y)=\sum_{x_i\in w} \sum_{y_i \in w}(f(x_i,y_i)-f(x_i-\Delta x,y_i-\Delta y))^2
C(x,y)=xi∈w∑yi∈w∑(f(xi,yi)−f(xi−Δx,yi−Δy))2
角点不会受光圈问题的影响，对于所有
Δ
x
,
Δ
y
，
S
w
(
Δ
x
,
Δ
y
)
\Delta x,\Delta y，S_w(\Delta x,\Delta y)
Δx,Δy，Sw(Δx,Δy)都是高响应，如果平移图像用一阶泰勒展开近似，则可表示为
f
(
x
i
−
Δ
x
,
y
i
−
Δ
y
)
≈
(
f
(
x
i
,
y
i
)
+
[
∂
f
(
x
i
,
y
i
)
∂
x
∂
f
(
x
i
,
y
i
)
∂
y
]
[
Δ
x
Δ
y
]
f(x_i-\Delta x,y_i-\Delta y)\approx (f(x_i,y_i) +\begin{bmatrix}\frac{\partial f(x_i,y_i)}{\partial x}&\frac{\partial f(x_i,y_i)}{\partial y}\end{bmatrix} \begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix}
f(xi−Δx,yi−Δy)≈(f(xi,yi)+[∂x∂f(xi,yi)∂y∂f(xi,yi)][ΔxΔy]
带入：
C
(
x
,
y
)
=
∑
x
i
∈
w
∑
y
i
∈
w
(
f
(
x
i
,
y
i
)
−
f
(
x
i
,
y
i
)
−
[
∂
f
(
x
i
,
y
i
)
∂
x
∂
f
(
x
i
,
y
i
)
∂
y
]
[
Δ
x
Δ
y
]
)
2
C(x,y)=\sum_{x_i\in w} \sum_{y_i \in w}\Bigg(f(x_i,y_i)-f(x_i,y_i)-\begin{bmatrix}\frac{\partial f(x_i,y_i)}{\partial x}&\frac{\partial f(x_i,y_i)}{\partial y}\end{bmatrix} \begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix}\Bigg)^2
C(x,y)=xi∈w∑yi∈w∑(f(xi,yi)−f(xi,yi)−[∂x∂f(xi,yi)∂y∂f(xi,yi)][ΔxΔy])2
=
∑
x
i
∈
w
∑
y
i
∈
w
(
−
[
∂
f
(
x
i
,
y
i
)
∂
x
∂
f
(
x
i
,
y
i
)
∂
y
]
[
Δ
x
Δ
y
]
)
2
=\sum_{x_i\in w} \sum_{y_i \in w}\Bigg(-\begin{bmatrix}\frac{\partial f(x_i,y_i)}{\partial x}&\frac{\partial f(x_i,y_i)}{\partial y}\end{bmatrix} \begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix}\Bigg)^2
=xi∈w∑yi∈w∑(−[∂x∂f(xi,yi)∂y∂f(xi,yi)][ΔxΔy])2
=
∑
x
i
∈
w
∑
y
i
∈
w
(
[
∂
f
(
x
i
,
y
i
)
∂
x
∂
f
(
x
i
,
y
i
)
∂
y
]
[
Δ
x
Δ
y
]
)
2
=\sum_{x_i\in w} \sum_{y_i \in w}\Bigg(\begin{bmatrix}\frac{\partial f(x_i,y_i)}{\partial x}&\frac{\partial f(x_i,y_i)}{\partial y}\end{bmatrix} \begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix}\Bigg)^2
=xi∈w∑yi∈w∑([∂x∂f(xi,yi)∂y∂f(xi,yi)][ΔxΔy])2
因为
u
2
=
u
T
u
u^2=u^Tu
u2=uTu
=
∑
x
i
∈
w
∑
y
i
∈
w
[
Δ
x
,
Δ
y
]
(
[
∂
f
∂
x
∂
f
∂
y
]
[
∂
f
∂
x
∂
f
∂
y
]
)
[
Δ
x
Δ
y
]
=\sum_{x_i\in w} \sum_{y_i \in w}[\Delta x,\Delta y]\Bigg(\begin{bmatrix} \frac {\partial f}{\partial x} \\ \frac {\partial f}{\partial y} \end{bmatrix}\begin{bmatrix} \frac {\partial f}{\partial x} \frac {\partial f}{\partial y} \end{bmatrix}\Bigg) \begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix}
=xi∈w∑yi∈w∑[Δx,Δy]([∂x∂f∂y∂f][∂x∂f∂y∂f])[ΔxΔy]
=
[
Δ
x
,
Δ
y
]
(
∑
x
i
∈
w
∑
y
i
∈
w
[
∂
f
∂
x
∂
f
∂
y
]
[
∂
f
∂
x
∂
f
∂
y
]
)
[
Δ
x
Δ
y
]
=[\Delta x,\Delta y]\Bigg(\sum_{x_i\in w} \sum_{y_i \in w}\begin{bmatrix} \frac {\partial f}{\partial x} \\ \frac {\partial f}{\partial y} \end{bmatrix}\begin{bmatrix} \frac {\partial f}{\partial x} \frac {\partial f}{\partial y} \end{bmatrix}\Bigg) \begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix}
=[Δx,Δy](xi∈w∑yi∈w∑[∂x∂f∂y∂f][∂x∂f∂y∂f])[ΔxΔy]
=
[
Δ
x
,
Δ
y
]
A
w
(
x
,
y
)
[
Δ
x
Δ
y
]
=[\Delta x,\Delta y]A_w(x,y)\begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix}
=[Δx,Δy]Aw(x,y)[ΔxΔy]
A
w
(
x
,
y
)
A_w(x,y)
Aw(x,y)表示图像区域
w
w
w在点
(
x
,
y
)
=
(
0
,
0
)
(x,y)=(0,0)
(x,y)=(0,0)处的二阶导数的一半。
A
A
A为：
A
(
x
,
y
)
=
2
∙
[
∑
x
i
∈
w
∑
y
i
∈
w
(
∂
f
(
x
i
,
y
i
)
∂
x
)
2
∑
x
i
∈
w
∑
y
i
∈
w
∂
f
(
x
i
,
y
i
)
∂
x
∂
f
(
x
i
,
y
i
)
∂
y
∑
x
i
∈
w
∑
y
i
∈
w
∂
f
(
x
i
,
y
i
)
∂
x
∂
f
(
x
i
,
y
i
)
∂
y
∑
x
i
∈
w
∑
y
i
∈
w
(
∂
f
(
x
i
,
y
i
)
∂
y
)
2
]
A(x,y)=2\bullet \begin{bmatrix} \sum_{x_i\in w} \sum_{y_i \in w}(\frac{\partial f(x_i,y_i)}{\partial x})^2 & \sum_{x_i\in w} \sum_{y_i \in w} \frac{\partial f(x_i,y_i)}{\partial x}\frac{\partial f(x_i,y_i)}{\partial y} \\ \sum_{x_i\in w} \sum_{y_i \in w} \frac{\partial f(x_i,y_i)}{\partial x}\frac{\partial f(x_i,y_i)}{\partial y}& \sum_{x_i\in w} \sum_{y_i \in w}(\frac{\partial f(x_i,y_i)}{\partial y})^2 \end{bmatrix}
A(x,y)=2∙[∑xi∈w∑yi∈w(∂x∂f(xi,yi))2∑xi∈w∑yi∈w∂x∂f(xi,yi)∂y∂f(xi,yi)∑xi∈w∑yi∈w∂x∂f(xi,yi)∂y∂f(xi,yi)∑xi∈w∑yi∈w(∂y∂f(xi,yi))2]
局部结构矩阵
A
A
A代表邻域的变化，其有两个特征值
λ
1
，
λ
2
\lambda_1，\lambda_2
λ1，λ2，有以下三种情形出现：
1、两个特征值都很小，这意味着，图像f在检测点处平坦——在该点没有边缘或角点；
2、一个特征值很小，另一个很大，局部邻域呈脊状，如有垂直于脊部的微小移动，图像f将发生显著的变化；
3、两个特征值都很大，任何方向的移动，都会造成图像f的显著变化，此时，检测到一个角点。
对于求解特征值需要一元二次方程，这必然会降低运算效率。因此，通过计算响应函数
R
(
A
)
=
d
e
t
(
A
)
−
k
(
t
r
a
c
e
(
A
)
)
2
R(A)=det(A)-k(trace(A))^2
R(A)=det(A)−k(trace(A))2，即
R
=
λ
1
λ
2
−
k
(
λ
1
+
λ
2
)
2
R=\lambda_1\lambda_2-k(\lambda_1+\lambda_2)^2
R=λ1λ2−k(λ1+λ2)2，
当
A
(
x
,
y
)
=
[
a
b
c
d
]
A(x,y)=\begin{bmatrix} a&b\\c&d \end{bmatrix}
A(x,y)=[acbd]
d
e
t
(
A
)
=
λ
1
λ
2
=
a
c
−
b
2
det(A)=\lambda_1\lambda_2=ac-b^2
det(A)=λ1λ2=ac−b2
t
r
a
c
e
(
A
)
=
λ
1
+
λ
2
=
a
+
c
trace(A)=\lambda_1+\lambda_2=a+c
trace(A)=λ1+λ2=a+c
检测步骤：
- 计算图像 I I I的梯度，获取水平和垂直方向上的梯度 I x I_x Ix和 I y I_y Iy；
- 构建协方差矩阵，计算每个像素点的协方差矩阵，即计算两个方向梯度的乘积 I x 2 = I x ∗ I x ， I y 2 = I y ∗ I y ， I x y = I x ∗ I y I_x^2=I_x*I_x，I_y^2=I_y*I_y，I_{xy}=I_x*I_y Ix2=Ix∗Ix，Iy2=Iy∗Iy，Ixy=Ix∗Iy；
- 使用加权窗口 w w w（可以是高斯加权，可以是平均加权），对 I x 2 , I y 2 , I x y I_x^2,I_y^2,I_{xy} Ix2,Iy2,Ixy进行加权处理，得到自相关元素 a = I x 2 ∗ w ， b = I x y 2 ∗ w ， c = I y 2 ∗ w a=I_x^2 *w，b=I_{xy}^2 *w，c=I_y^2 *w a=Ix2∗w，b=Ixy2∗w，c=Iy2∗w；
- 计算协方差矩阵的特征值 R R R，判断每个像素点是否为角点。
- 对角点进行非极大值抑制，去除重复的角点。
- 根据设定的阈值，筛选出最终的角点。
OpenCV函数
void cv::cornerHarris(InputArray src, OutputArray dst, int blockSize, int ksize, double k, int borderType=BORDER_DEFAULT)
Parameters
src 输入图像，单通道8位或浮点型；
dst 输出Harris算子响应图像，该变量存储的R值，类型为单通道32位浮点型；
blockSize 步骤3中的加权窗口尺寸；
ksize Sobel算子求梯度的尺寸；
k 响应函数中的k；
borderType 填充边缘方法，见边缘填充博客；
OpenCV源码分析
等待后续…
