Harris角点检测器
1 知识点
角点定义:
a. 局部窗口沿各方向移动,均产生明显变化的点
b. 图像局部曲率突变的点
角点图释:
常见的角点:
相关名词:
两个像素块(大小相同):
I
1
(
x
)
I_1(x)
I1(x),
I
2
(
x
)
I_2(x)
I2(x)
窗口函数:
Σ
x
,
y
w
(
x
,
y
)
\Sigma_{x,y}w(x,y)
Σx,yw(x,y)
灰度变化:
E
(
u
,
v
)
E(u,v)
E(u,v)
角点响应函数:
R
R
R
权重矩阵:
W
W
W(通常为高斯滤波器
G
σ
G_\sigma
Gσ)
相关计算:
-
将窗口平移 [ u , v ] [u,v] [u,v]产生灰度变化 E ( u , v ) E(u,v) E(u,v)
E ( u , v ) = Σ x , y w ( x , y ) [ I ( x + u , y + v ) − I ( x , y ) ] E(u,v)=\Sigma_{x,y}w(x,y)[I(x+u,y+v)-I(x,y)] E(u,v)=Σx,yw(x,y)[I(x+u,y+v)−I(x,y)]
利用一级泰勒展开式,可得,
I ( x + u , y + v ) = I ( x , y ) + I x u + I y v + O ( u 2 , v 2 ) I(x+u,y+v)=I(x,y)+I_xu+I_yv+O(u^2,v^2) I(x+u,y+v)=I(x,y)+Ixu+Iyv+O(u2,v2)
E ( u , v ) ≈ Σ x , y w ( x , y ) [ I x u + I y v ] 2 ≅ [ u v ] M [ u v ] E(u,v)≈\Sigma_{x,y}w(x,y)[I_xu+I_yv]^2\cong \begin{bmatrix}u \ v\end{bmatrix}M \begin{bmatrix}u \\ v\end{bmatrix} E(u,v)≈Σx,yw(x,y)[Ixu+Iyv]2≅[u v]M[uv]
其中,M = Σ x , y w ( x , y ) [ I x 2 I x I y I x I y I y 2 ] M=\Sigma_{x,y}w(x,y)\begin{bmatrix} I_x^2 & I_xI_y \\ I_xI_y & I_y^2 \\ \end{bmatrix} M=Σx,yw(x,y)[Ix2IxIyIxIyIy2]
-
点 x x x的对称半正定矩阵:
M I = M I ( x ) = ∇ I ∇ I T = [ I x I y ] [ I x I y ] = [ I x 2 I x I y I x I y I y 2 ] M_I=M_I(x)=\nabla I \nabla I^T=\begin{bmatrix} I_x \\ I_y \\ \end{bmatrix} \begin{bmatrix} I_x \ I_y\\ \end{bmatrix} =\begin{bmatrix} I_x^2 & I_xI_y \\ I_xI_y & I_y^2 \\ \end{bmatrix} MI=MI(x)=∇I∇IT=[IxIy][Ix Iy]=[Ix2IxIyIxIyIy2]
M I M_I MI周围像素上的局部平均为:
M I ˉ = W ∗ M I \bar{M_I}=W*M_I MIˉ=W∗MI
通过式 λ E − M I = 0 \lambda E-M_I=0 λE−MI=0,得到特征值 λ 1 = ∣ ∇ I ∣ 2 \lambda_1=|\nabla I|^2 λ1=∣∇I∣2, λ 2 = 0 \lambda_2=0 λ2=0
特征值会依赖于局部图像变化,具体表现为:图像梯度变化 ⟹ \implies ⟹ M I ˉ \bar{M_I} MIˉ特征值变化。
不同的特征值对应不同的像素情况:
a. λ 1 \lambda_1 λ1、 λ 2 \lambda_2 λ2都是很大的正数 ⟹ \implies ⟹ x x x为角点
b. λ 1 \lambda_1 λ1很大, λ 2 = 0 \lambda_2=0 λ2=0 ⟹ \implies ⟹区域存在一个边, M I ˉ \bar{M_I} MIˉ变化不大
c. λ 1 = λ 2 = 0 \lambda_1= \lambda_2=0 λ1=λ2=0 ⟹ \implies ⟹区域为空
- 角点响应
R = d e t M I − k ( t r a c e M I ) 2 R = det M_I-k(trace M_I)^2 R=detMI−k(traceMI)2或
R = d e t M I ( t r a c e M I ) 2 R = \frac{det M_I}{(trace M_I)^2} R=(traceMI)2detMI
其中,
d e t M I = λ 1 λ 2 det M_I = \lambda_1 \lambda_2 detMI=λ1λ2
t r a c e M I = λ 1 + λ 2 trace M_I = \lambda_1+\lambda_2 traceMI=λ1+λ2
R 只与M的特征值有关
角点:R 为大数值正数
边缘:R为大数值负数
平坦区:R为小数值
-
角点进行挑选方法:
a. 设置阈值。选取像素值高于阈值的所有图像点
b. 角点之间间隔必须大于设定的最小距离 -
在图像中寻找对应点:
Harris角点检测器仅能够检测出图像中的兴趣点,但没有给出通过比较图像间的兴趣点来寻找匹配角点的方法。
Harris角点的描述子通常是由周围图像像素块的灰度值以及用于比较的归一化互相关矩阵构成的。
两个像素块 I 1 ( x ) I_1(x) I1(x), I 2 ( x ) I_2(x) I2(x)的相关矩阵为
c ( I 1 , I 2 ) = Σ x f ( I 1 ( x ) , I 2 ( x ) ) , f 指 的 是 相 关 方 法 c(I_1,I_2)=\Sigma_{x}f(I_1(x),I_2(x)),f指的是相关方法 c(I1,I2)=Σxf(I1(x),I2(x)),f指的是相关方法
互相关矩阵为
f ( I 1 , I 2 ) = I 1 I 2 f(I_1,I_2)=I_1I_2 f(I1,I2)=I1I2
∴ c ( I 1 , I 2 ) = I 1 ⋅ I 2 \therefore c(I_1,I_2)=I_1·I_2 ∴c(I1,I2)=I1⋅I2
⋅ · ⋅ 指向量乘法, c ( I 1 , I 2 ) c(I_1,I_2) c(I1,I2)值越大, I 1 , I 2 I_1,I_2 I1,I2相似度越高。
归一化互相关矩阵为:
n c c ( I 1 , I 2 ) = 1 n − 1 Σ x ( I 1 ( x ) − μ 1 ) σ 1 ⋅ ( I 2 ( x ) − μ 2 ) σ 2 ncc(I_1,I_2)=\frac{1}{n-1}\Sigma_{x}\frac{(I_1(x)-\mu_1)}{\sigma_1}·\frac{(I_2(x)-\mu_2)}{\sigma_2} ncc(I1,I2)=n−11Σxσ1(I1(x)−μ1)⋅σ2(I2(x)−μ2)
n n n指像素块中像素的数目,像素块即以该像素点为中心的周围矩形部分图像
μ \mu μ指平均像素强度
σ \sigma σ指像素块的标准差
2 实验部分
- 角点检测
a. 返回角点响应值 R R R
def compute_harris_response(im,sigma=3):
""" Compute the Harris corner detector response function
for each pixel in a graylevel image. """
# derivatives
imx = zeros(im.shape)
filters.gaussian_filter(im, (sigma,sigma), (0,1), imx)
imy = zeros(im.shape)
filters.gaussian_filter(im, (sigma,sigma), (1,0), imy)
# compute components of the Harris matrix
Wxx = filters.gaussian_filter(imx*imx,sigma)
Wxy = filters.gaussian_filter(imx*imy,sigma)
Wyy = filters.gaussian_filter(imy*imy,sigma)
# determinant and trace
Wdet = Wxx*Wyy - Wxy**2
Wtr = Wxx + Wyy
return Wdet / Wtr
b. 返回角点并标记
def get_harris_points(harrisim,min_dist=10,threshold=0.1):
""" Return corners from a Harris response image
min_dist is the minimum number of pixels separating
corners and image boundary. """
# find top corner candidates above a threshold
corner_threshold = harrisim.max() * threshold
harrisim_t = (harrisim > corner_threshold) * 1
# get coordinates of candidates
coords = array(harrisim_t.nonzero()).T
# ...and their values
candidate_values = [harrisim[c[0],c[1]] for c in coords]
# sort candidates (reverse to get descending order)
index = argsort(candidate_values)[::-1]
# store allowed point locations in array
allowed_locations = zeros(harrisim.shape)
allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1
# select the best points taking min_distance into account
filtered_coords = []
for i in index:
if allowed_locations[coords[i,0],coords[i,1]] == 1:
filtered_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 filtered_coords
def plot_harris_points(image,filtered_coords):
""" Plots corners found in image. """
figure()
gray()
imshow(image)
plot([p[1] for p in filtered_coords],
[p[0] for p in filtered_coords],'ro')
axis('off')
show()
c. 运行主函数
# 读入图像
im = array(Image.open(r"C:\Users\13121\Desktop\test.jpg").convert('L'))
# 检测harris角点
harrisim = compute_harris_response(im)
# Harris响应函数
harrisim1 = 255 - harrisim
figure()
gray()
#画出Harris响应图
subplot(141)
imshow(harrisim1)
axis('off')
axis('equal')
threshold = [0.01, 0.05, 0.1]
for i, thres in enumerate(threshold):
filtered_coords = get_harris_points(harrisim, 6, thres)
subplot(1, 4, i+2)
imshow(im)
plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '.')
axis('off')
show()
d. 运行结果图
域值越大,得到的角点越稀疏
也可以直接调用函数得到角点,结果图为:
代码为:
im = array(Image.open(r"C:\Users\13121\Desktop\test.jpg").convert('L'))
harrisim = compute_harris_response(im)
filter_coords = get_harris_points(harrisim,6)
plot_harris_points(im, filtered_coords)
- 图像中找对应点
a. 选用奇数像素的灰度图像块
def get_descriptors(image,filtered_coords,wid=5):
""" For each point return pixel values around the point
using a neighbourhood of width 2*wid+1. (Assume points are
extracted with min_distance > wid). """
desc = []
for coords in filtered_coords:
patch = image[coords[0]-wid:coords[0]+wid+1,
coords[1]-wid:coords[1]+wid+1].flatten()
desc.append(patch)
return desc
b. 归一化互相关
def match(desc1,desc2,threshold=0.5):
""" For each corner point descriptor in the first image,
select its match to second image using
normalized cross correlation. """
n = len(desc1[0])
# pair-wise distances
d = -ones((len(desc1),len(desc2)))
for i in range(len(desc1)):
for j in range(len(desc2)):
d1 = (desc1[i] - mean(desc1[i])) / std(desc1[i])
d2 = (desc2[j] - mean(desc2[j])) / std(desc2[j])
ncc_value = sum(d1 * d2) / (n-1)
if ncc_value > threshold:
d[i,j] = ncc_value
ndx = argsort(-d)
matchscores = ndx[:,0]
return matchscores
c. 两边对称匹配
def match_twosided(desc1,desc2,threshold=0.5):
""" Two-sided symmetric version of match(). """
matches_12 = match(desc1,desc2,threshold)
matches_21 = match(desc2,desc1,threshold)
ndx_12 = where(matches_12 >= 0)[0]
# remove matches that are not symmetric
for n in ndx_12:
if matches_21[matches_12[n]] != n:
matches_12[n] = -1
return matches_12
d. 两幅图像拼接
def appendimages(im1,im2):
""" Return a new image that appends the two images side-by-side. """
# select the image with the fewest rows and fill in enough empty rows
rows1 = im1.shape[0]
rows2 = im2.shape[0]
if rows1 < rows2:
im1 = concatenate((im1,zeros((rows2-rows1,im1.shape[1]))),axis=0)
elif rows1 > rows2:
im2 = concatenate((im2,zeros((rows1-rows2,im2.shape[1]))),axis=0)
# if none of these cases they are equal, no filling needed.
return concatenate((im1,im2), axis=1)
def plot_matches(im1,im2,locs1,locs2,matchscores,show_below=True):
""" Show a figure with lines joining the accepted matches
input: im1,im2 (images as arrays), locs1,locs2 (feature locations),
matchscores (as output from 'match()'),
show_below (if images should be shown below matches). """
im3 = appendimages(im1,im2)
if show_below:
im3 = vstack((im3,im3))
imshow(im3)
cols1 = im1.shape[1]
for i,m in enumerate(matchscores):
if m>0:
plot([locs1[i][1],locs2[m][1]+cols1],[locs1[i][0],locs2[m][0]],'c')
axis('off')
e. 运行结果:
im1 = array(Image.open(r"C:\Users\13121\Desktop\test.jpg").convert("L"))
im2 = array(Image.open(r"C:\Users\13121\Desktop\test.jpg").convert("L"))
wid = 5
harrisim = compute_harris_response(im1, 5)
filtered_coords1 = get_harris_points(harrisim, wid+1)
d1 = get_descriptors(im1, filtered_coords1, wid)
harrisim = compute_harris_response(im2, 5)
filtered_coords2 = get_harris_points(harrisim, wid+1)
d2 = get_descriptors(im2, filtered_coords2, wid)
print('starting matching')
matches = match_twosided(d1, d2)
figure()
gray()
plot_matches(im1, im2, filtered_coords1, filtered_coords2, matches)
show()
图形像素块的互相关矩阵有较弱的描述性,Harris不具有尺度不变性和旋转不变性,算法中的像素块的大小也会影响对应匹配的结果。