阿齐兹的PyCV学习笔记——Harris角点检测器

Harris角点检测器

Harris角点检测算法：如果像素周围显示存在多于一个方向的边，就认为该点是兴趣点，该点就称为角点

对图像中的每一个像素，计算Harris矩阵，Harris矩阵的特征值有三种情况：

都是很大的正数，则该点是角点

一个很大，一个约等于0，则该区域存在一个边

都约等于0，该区域为空

在不需要计算特征值的情况下，引入指示函数区分重要的情况

定义角点响应函数，参数sigma定义了高斯滤波器的尺度大小，返回像素值为Harris响应函数的一幅图像。需要使用scipy.ndimage.filters模块中的高斯导数滤波器计算导数，因为在角点检测中需要抑制噪声

from scipy.ndimage import filters
def compute_harris_response(im,sigma=3):
"""在一幅灰度图像中，对每个像素计算Harris角点检测器响应函数"""

#计算导数
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)

#计算Harris矩阵分量
Wxx = filters.gaussian_filter(imx*imx,sigma)
Wxy = filters.gaussian_filter(imx*imy,sigma)
Wyy = filters.gaussian_filter(imy*imy,sigma)

#计算特征值和迹
Wdet = Wxx*Wyy - Wxy**2
Wtr = Wxx + Wyy

return Wdet / Wtr

定义函数get_harris_points()，选取像素值高于阈值的所有图像点，再加上额外的限制：角点之间的间隔必须大于设定的最小距离，返回以角点响应值递减的顺序排序的角点

def get_harris_points(harrisim,min_dist=10,threshold=0.1):
"""从一幅Harris响应图像中返回角点，min_dist为分割角点和图像边界的最少像素数目"""

#寻找高于阈值的候选角点
corner_threshold = harrisim.max()*threshold
harrisim_t = (harrisim>corner_threshold)*1

#得到候选点的坐标
coords = array(harrisim_t.nonzero()).T

#以及他们的Harris响应值
candidate_values = [harrisim[c[0],c[1]] for c in coords]

#对候选点按照Harris响应值进行排序
index = argsort(candidate_values)

#将可行点的位置保存到数组中
allowed_locations = zeros(harrisim.shape)
allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1

#按照min_distance原则，选择最佳Harris点
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

使用Matplotlib模块绘制角点

def plot_harris_points(image,filtered_coords):
"""绘制图像中检测到的角点"""

figure()
gray()
imshow(image)
plot([p[1] for p in filtered_coords],[p[0] for p in filtered_coords],'*')
axis('off')
show()

检测角点的测试代码：

im = array(Image.open('C:/Users/0AQZ0/Documents/exercisecode/Python/PyCV/Images/001.jpg').convert('L'))
harrisim = compute_harris_response(im)
filtered_coords = get_harris_points(harrisim,6)
plot_harris_points(im,filtered_coords)

效果图：



兴趣点描述子是分配给兴趣点的一个向量，描述该点附近图像的表观信息。描述子越好，寻找到的对应点越好。对应点是指描述相同物体和场景在不同图像上形成的像素点

Harris角点的描述子通常是周围图像像素块的灰度值，以及用于比较的归一化互相关矩阵构成的。图像的像素块由以该像素点为中心的周围矩阵部分图像构成

get_descriptors()方法可以将图像块像素值压平成一个向量，然后添加到描述子列表中，第一个参数是奇数大小长度的方形灰度图像块，图像块的中心是处理的像素点

def get_descriptors(image,filtered_coords,wid=5):
"""对每个返回的点，返回点周围2*wid+1个像素的值（假设选取点的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

match()方法使用归一化的互相关矩阵，将每个描述子匹配到另一个图像中的最优的候选点

def match(desc1,desc2,threshold=0.5):
"""对于第一幅图像中的每个角点描述子，使用归一化互相关，选取它在第二幅图像中的匹配焦点"""

n=len(desc1[0])

#点对的距离
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

match_twosided()方法可以从第二幅图像向第一幅图像匹配，过滤掉在两种方法中不都是最好的匹配

def match_twosided(desc1,desc2,threshold=0.5):
"""两边对称版本的match()"""

matches_12 = match(desc1,desc2,threshold)
matches_21 = match(desc2,desc1,threshold)

ndx_12 = where(matches_12 >= 0)[0]

#去除非对称的匹配
for n in ndx_12:
if matches_21[matches_12
] != n:
matches_12
= -1
return matches_12

appendimages()方法、plot_matches()方法可以实现匹配点的可视化，绘制图像使用线段连接匹配的像素点

def appendimages(im1,im2):
"""返回将两幅图像并排拼接成的一幅新图像"""

#选取具有最少行数的图像，然后填充足够的空行
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)
#如果这些情况都没有，那么它们的行数相同，不需要进行填充

return concatenate((im1,im2),axis=1)

def plot_matches(im1,im2,locs1,locs2,matchscores,show_below=True):
"""显示一幅带有连接匹配之间连线的图片
输入：im1，im2（数组图像），locs1，locs2（特征位置），matchscores（match（）的输出），
show_below（如果图像应该显示在匹配的下方）"""

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')

使用归一化互相关矩阵寻找对应点的实例：

im1 = array(Image.open('C:/Users/0AQZ0/Desktop/中控杯/视觉识别代码/images/01.jpg').convert('L'))
im2 = array(Image.open('C:/Users/0AQZ0/Desktop/中控杯/视觉识别代码/images/02.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 'start matching'
matches = match_twosided(d1,d2)

figure()
gray()
plot_matches(im1,im2,filtered_coords1,filtered_coords2,matches)
show()

效果图：



缺点：算法的结果存在一些不正确匹配，因为图像像素块的互相关矩阵具有较弱的描述性，而且描述符不具有尺度不变性和旋转不变性，算法中像素块的大小也会影响匹配的结果