Harris角点检测
Harris角点检测
角点的定义:
角点检测算法基本思想是使用一个固定窗口(取某个像素的一个邻域窗口)在图像上进行任意方向上的滑动,比较滑动前与滑动后两种情况,窗口中的像素灰度变化程度,如果存在任意方向上的滑动,都有着较大灰度变化,那么我们可以认为该窗口中存在角点。
角点:最直观的印象就是在水平和竖直两个方向变化均较大的两个点,即X,Y方向梯度都较大
边缘:仅在水平或竖直方向有较大的变化量,即X、Y方向梯度只有一个较大
平坦区域:在水平和竖直方向的变化量均较小,即X、Y方向梯度都较小
在特征点检测中经常提出尺度不变、旋转不变、抗噪声影响等,这些是判断特征点是否稳定的指标。
性能较好的角点:
- 检测出图像中“真实”的角点
- 准确的定位性能
- 很高的重复检测率
- 噪声的鲁棒性
- 较高的计算效率
今天看到角点检测的时候真的是一脸蒙蔽,服了,书上直接给出的公式。整的我直接懵了。参考晚上其他大佬写的博客后逐渐明白了一点
书上给的公式如下:
E
(
u
,
v
)
=
∑
x
,
y
w
(
x
,
y
)
[
]
I
x
2
I
x
I
y
I
x
I
y
I
y
2
]
(B)
E(u,v)=\sum_{x,y}w(x,y) \left[] \begin{matrix} Ix^2 & IxIy\\ IxIy & Iy^2 \end{matrix} \right] \tag{B}
E(u,v)=x,y∑w(x,y)[]Ix2IxIyIxIyIy2](B)
w(x,y)是窗口函数,最简单情形就是窗口W内的所有像素所对应的w权重系数均为1.但有时候,我们会将w(x,y)函数设置为以窗口W中心为原点的二元正态分布。如果窗口W中心点是角点时,移动前与移动后,该点在灰度变化贡献最大;而离窗口W中心(角点)较远的点,这些点的灰度变化几近平缓,这些点的权重系数,可以设定小值,以示该点对灰度变化贡献较小,那么我们自然而然想到使用二元高斯函数来表示窗口函数;
Harris矩阵推导过程:
使用一个固定窗口在图像上进行任意方向上的滑动。滑动分别为u和v.
比较滑动前与滑动后俩种情况。所以可以得到以下公式:
E
(
u
,
v
)
=
∑
x
,
y
w
(
x
,
y
)
[
I
(
x
+
u
,
y
+
v
)
−
I
(
x
,
y
)
]
2
E(u,v)=\sum_{x,y}w(x,y)[I(x+u,y+v)-I(x,y)]^2
E(u,v)=x,y∑w(x,y)[I(x+u,y+v)−I(x,y)]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)+ufx(x,y)+vfy(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
=
∑
[
I
(
x
,
y
)
+
u
I
x
+
v
I
y
−
I
(
x
,
y
)
]
2
所以\sum[I(x+u,y+v)-I(x,y)]^2=\sum[I(x,y)+uIx+vIy-I(x,y)]^2
所以∑[I(x+u,y+v)−I(x,y)]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
=
∑
[
u
,
v
]
[
I
x
2
I
x
I
y
I
x
I
y
I
y
2
]
[
u
v
]
(B)
=\\\sum(u^2Ix^2)+2uvIxIy+v^2Iy^2=\\\sum[u,v] \left[ \begin{matrix} Ix^2 & IxIy\\ IxIy & Iy^2 \end{matrix} \right] \tag{B} \left[ \begin{matrix} u\\ v\\ \end{matrix} \right]
=∑(u2Ix2)+2uvIxIy+v2Iy2=∑[u,v][Ix2IxIyIxIyIy2][uv](B)
- [u,v]是窗口W的偏移量;
- (x,y)是窗口W所对应的像素坐标位置,窗口有多大,就有多少个位置;
- I(x,y)是像素坐标位置(x,y)的图像灰度值;
- I(x+u,y+v)是像素坐标位置(x+u,y+v)的图像灰度值;
Ix,Iy分别为窗口内像素点(x,y)在x方向上和y方向上的梯度值。 其中**窗口函数(权重矩阵)**可以是平坦的,通常为高斯滤波器Gσ:
设:
M
i
=
[
I
x
2
I
x
I
y
I
x
I
y
I
y
2
]
M
=
∑
x
,
y
w
(
x
,
y
)
∗
M
i
Mi =\left[ \begin{matrix} Ix^2 & IxIy\\ IxIy & Iy^2 \end{matrix} \right]\\ M = \sum_{x,y}w(x,y) * Mi
Mi=[Ix2IxIyIxIyIy2]M=x,y∑w(x,y)∗Mi
所以可以得到卷积:M
该卷积的目的是得到MI在周围像素上的局部平均。**矩阵M又称为Harris矩阵。**W 的宽度决定了在像素x 周围的感兴趣区域。
根据上述表达式,当窗口在平坦区域上移动,可以想象得到,灰度不会发生什么变换。E(u,v)=0;如果窗口处在纹理比较丰富的区域上滑动,那么灰度变化会很大。算法最终思想就是计算灰度发生较大变化时所对应的位置,当然这个较大是指任意方向上的滑动,并非单指某个方向。
接下来我们对矩阵M的特征值进行分析。设矩阵M的特征值为λ1,λ2.
则有以下三种情况:
具体为什么有这三种情况以下就是推导过程:
我们可以将E(u,v)近似为二项函数:
E
(
u
,
v
)
=
A
u
2
+
2
C
u
v
+
B
v
2
=
[
A
C
C
A
]
A
=
∑
x
,
y
w
(
x
,
y
)
∗
I
x
2
B
=
∑
x
,
y
w
(
x
,
y
)
∗
I
y
2
C
=
∑
x
,
y
w
(
x
,
y
)
∗
I
x
I
y
E(u,v)=Au^2+2Cuv+Bv^2=\left[ \begin{matrix} A & C\\ C & A\\ \end{matrix} \right]\\ A=\sum_{x,y}w(x,y)*Ix^2\\ B=\sum_{x,y}w(x,y)*Iy^2\\ C=\sum_{x,y}w(x,y)*IxIy
E(u,v)=Au2+2Cuv+Bv2=[ACCA]A=x,y∑w(x,y)∗Ix2B=x,y∑w(x,y)∗Iy2C=x,y∑w(x,y)∗IxIy
二次项函数本质上就是一个椭圆函数。椭圆的长和宽是由MM的特征值λ1,λ2决定的(椭圆的长短轴正是矩阵MM特征值平方根的倒数),椭圆的方向是由M的特征向量决定的,椭圆方程为:
[
u
v
]
∗
M
∗
[
u
v
]
=
1
\left[ \begin{matrix} u & v\\ \end{matrix} \right]*M*\left[ \begin{matrix} u\\ v\\ \end{matrix} \right]=1
[uv]∗M∗[uv]=1
虽然我们利用E(u,v)来描述角点的基本思想,然而最终我们仅仅使用的是矩阵M。让我们看看矩阵M形式,是不是跟协方差矩阵形式很像,像归像,但是还是有些不同,哪儿不同?一般协方差矩阵对应维的随机变量需要减去该维随机变量的均值:
把Ix,Iy看成两个字段,假设窗口内有m个像素点,也就是等价于有m个样本,我们先计算每个字段的均值:
I
x
‾
=
∑
i
m
I
x
i
I
y
‾
=
∑
i
m
I
y
i
\overline{Ix}=\sum_{i}^mIxi\\ \overline{Iy}=\sum_{i}^mIyi
Ix=i∑mIxiIy=i∑mIyi
我们仍然使用(Ixi,Iyi)表示样本(Ixi,Iyi)去均值后的值,则由这m个样本组成的矩阵:
x
=
[
I
x
1
I
x
2
.
.
.
I
x
m
I
y
1
I
y
2
.
.
.
I
y
m
]
x=\left[ \begin{matrix} Ix1 & Ix2 & ... & Ixm \\ Iy1 & Iy2 & ... & Iym\\ \end{matrix} \right]
x=[Ix1Iy1Ix2Iy2......IxmIym]
则对应协方差矩阵可以写成:
C
=
1
m
∗
X
∗
X
T
=
1
m
∗
∑
i
=
1
[
I
x
2
I
x
I
y
I
x
I
y
I
y
2
]
C=\frac{1}{m}*X*X^T=\frac{1}{m}*\sum_{i=1}\left[ \begin{matrix} Ix^2 & IxIy\\ IxIy & Iy^2 \end{matrix} \right]
C=m1∗X∗XT=m1∗i=1∑[Ix2IxIyIxIyIy2]
但矩阵M中并没有这样做,所以在矩阵M里,我们先进行各维的均值化处理,那么各维所对应的随机变量的均值为0,协方差矩阵就大大简化了,简化的最终结果就是矩阵M,是否明白了(注意为了简化运算,我们先假设M矩阵中的权重系数w(x,y)=1,并且省略掉了求均值.我们的目的是分析数据的主要成分
如果我们对协方差矩阵M进行对角化,很明显,其对角线就是各个字段的方差,这点大家应该明白吧?不明白的话可以复习下PCA原理。
- 如果两个字段(Ix,Iy)所对应的特征值都比较大,说明什么? 像素点的梯度分布比较散,梯度变化程度比较大,符合角点在窗口区域的特点;
- 如果是平坦区域,那么像素点的梯度所构成的点集比较集中在原点附近,因为窗口区域内的像素点的梯度幅值非常小,此时矩阵M的对角化的两个特征值比较小;
- 如果是边缘区域,在计算像素点的x、y方向上的梯度时,边缘上的像素点的某个方向的梯度幅值变化比较明显,另一个方向上的梯度幅值变化较弱,其余部分的点都还是集中原点附近,这样MM对角化后的两个特征值理论应该是一个比较大,一个比较小,当然对于边缘这种情况,可能是呈45°的边缘,致使计算出的特征值并不是都特别的大,总之跟含有角点的窗口的分布情况还是不同的。
因此可以得出下列结论:
-
特征值都比较大时,即窗口中含有角点;
-
特征值一个较大,一个较小,窗口中含有边缘;
-
特征值都比较小,窗口处在平坦区域;
那么根据书上又提出了一个不用实际计算特征值的方法。也就是度量角点响应:
对每一个窗口计算得到一个分数R,根据R的大小来判定窗口内是否存在harris特征角。分数R根据下面公式计算得到:
R = d e t ( M ) − k ( t r a c e ( M ) ) 2 其 中 : d e t ( M ) = λ 1 ∗ λ 2 t r a c e ( M ) = λ 1 + λ 2 另 外 我 在 别 的 地 方 看 到 简 化 后 的 R 为 : R = [ ( I x 2 + I y 2 ) − ( I x I y ) 2 ] / ( I x 2 + I y 2 ) R=det(M)-k(trace(M))^2\\ 其中:det(M)=\lambda1*\lambda2\\ trace(M)=\lambda1+\lambda2\\ 另外我在别的地方看到简化后的R为: R=[(Ix^2+Iy^2)-(IxIy)^2]/(Ix^2+Iy^2) R=det(M)−k(trace(M))2其中:det(M)=λ1∗λ2trace(M)=λ1+λ2另外我在别的地方看到简化后的R为:R=[(Ix2+Iy2)−(IxIy)2]/(Ix2+Iy2)
得到Harris角点的步骤
所以得到Harris角点的步骤为:
1、计算图像I(x,y)I(x,y)在x和y两个方向的梯度Ix,Iy;
2、计算图像两个方向梯度的乘积;
3、使用高斯函数对第二步算出的数进行高斯加权,计算中心点为(x,y)的窗口W对应的矩阵M
4、计算每个像素点(x,y)处的(x,y)处的Harris响应值R;
5:过滤大于某一阈值t的R值;
代码实现
python代码如下:
from scipy.ndimage import filters
from PIL import Image
from pylab import *
from imtools import imresize
np.set_printoptions(threshold=np.inf)#解决printf打印矩阵不完全的问题
def harris_response(im,sigma=3):
"""计算图像的harris响应函数"""
#计算导数
imx = zeros(im.shape)
imx = filters.gaussian_filter(im,(sigma,sigma),(0,1),imx)
imy = zeros(im.shape)
imy = filters.gaussian_filter(im, (sigma, sigma), (1, 0), imy)
#计算Harris的各个分量
wxx = filters.gaussian_filter(imx*imx,sigma)
wxy = filters.gaussian_filter(imx*imy,sigma)
wyy = filters.gaussian_filter(imy*imy,sigma)
#计算像素的角点响应函数
return (wxx*wyy - 2*wxy)/(wxx + wyy)
def get_harris_points(harrism,min_dist = 10,thresold = 0.1):
"""从一幅Harrisim响应中返回角点,min_dist为分割角点和图像边界的最少像素数目"""
corner_thsold = harrism.max()*thresold
harrism_t = (harrism > corner_thsold) * 1
#得到候选点的坐标
coords = array(harrism_t.nonzero()).T#返回非零值的坐标的矩阵
#他们的Harris响应值
candidate_values = [harrism[c[0],c[1]] for c in coords]
#对候选点进行harris响应值进行排序
index = argsort(candidate_values)[::-1]#将x中的元素从小到大排列,提取其对应的index(索引),然后输出到y
#将可行点的位置保存在数组里
allowed_locations = zeros(harrism.shape)
allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1
#按照min_distance原则,选择最佳harris点
filters_coords = []
for i in index:
if allowed_locations[coords[i,0],coords[i,1]] == 1:
filters_coords.append(coords[i])
allowed_locations[(coords[i,0]-min_dist):(coords[i,0]+min_dist),(coords[i,1]-min_dist):coords[i,1]+min_dist] = 0
return filters_coords
本文来自博客园,作者:Bathwind_W,转载请注明原文链接:https://www.cnblogs.com/bathwind/p/18107960