文章目录

1.Harris角点检测1.1基础知识1.2算法思想1.3Harris数学公式1.4代码实例效果

1.Harris角点检测

1.1基础知识

使用一个滑动窗口在下面三幅图中滑动，可以得出以下结论：

左图表示一个平坦区域，在各方向移动，窗口内像素值均没有太大变化；

中图表示一个边缘特征（Edges），如果沿着水平方向移动(梯度方向)，像素值会发生跳变；如果沿着边缘移动(平行于边缘) ，像素值不会发生变化；

右图表示一个角（Corners），不管你把它朝哪个方向移动，像素值都会发生很大变化。

1.2算法思想

算法的核心是利用局部窗口在图像上进行移动，判断灰度是否发生较大的变化。如果窗口内的灰度值（在梯度图上）都有较大的变化，那么这个窗口所在区域就存在角点。

这样就可以将 Harris 角点检测算法分为以下三步：

当窗口（局部区域）同时向 x （水平）和 y（垂直） 两个方向移动时，计算窗口内部的像素值变化量 E(x,y) ；对于每个窗口，都计算其对应的一个角点响应函数 R；然后对该函数进行阈值处理，如果 R>threshold，表示该窗口对应一个角点特征。

1.3Harris数学公式

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

[u,v][u,v]是窗口WW的偏移量；(x,y)(x,y)是窗口WW所对应的像素坐标位置，窗口有多大，就有多少个位置；I(x,y)是像素坐标位置(x,y)的图像灰度值；I(x+u,y+v)是像素坐标位置(x+u,y+v)的图像灰度值；

1.4代码实例

from PIL import Imagefrom pylab import *from scipy.ndimage import filtersdef compute_harris_response(img, sigma=3):# 计算导数imx = zeros(img.shape)filters.gaussian_filter(img, (sigma, sigma), (0, 1), imx)imy = zeros(img.shape)filters.gaussian_filter(img, (sigma, sigma), (1, 0), imy)# 可视化原图plt.subplot(2, 3, 1)plt.title('灰度图')plt.axis('off')plt.imshow(img, plt.cm.gray)# 可视化x方向导数图像plt.subplot(2, 3, 2)plt.title('x方向导数')plt.axis('off')plt.imshow(imx, plt.cm.gray)# 可视化y方向导数图像plt.subplot(2, 3, 3)plt.title('y方向导数')plt.axis('off')plt.imshow(imy, plt.cm.gray)plt.subplot(2, 3, 4)plt.title('x方向二阶导数')plt.axis('off')plt.imshow(imx * imx, plt.cm.gray)plt.subplot(2, 3, 5)plt.title('x/y方向二阶导数')plt.axis('off')plt.imshow(imx * imy, plt.cm.gray)plt.subplot(2, 3, 6)plt.title('y方向二阶导数')plt.axis('off')plt.imshow(imy * imy, plt.cm.gray)plt.show()# 计算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 + Wyyreturn Wdet / Wtr# 从一幅Harris响应图像中返回角点。min_dist为分割角点和图像边界的最少像素数目def get_harris_points(harrisim, min_dist=10, threshold=0.1):# 寻找高于阈值的候选角点corner_threshold = harrisim.max() * threshold# print(corner_threshold)harrisim_t = (harrisim > corner_threshold) * 1# plt.imshow(harrisim_t, plt.cm.gray)# 得到候选点的坐标(取出是角点的坐标)coords = array(harrisim_t.nonzero()).T# print(coords.shape)# 以及它们的Harris响应值candidate_values = [harrisim[c[0], c[1]] for c in coords]# 对候选点按照Harris响应值进行排序index = argsort(candidate_values)# 将可行点的位置保存到数组中allowed_location = zeros(harrisim.shape)# 相当于把外面一圈有响应值的点先舍去（匹配主体在中心）allowed_location[min_dist:-min_dist, min_dist:-min_dist] = 1# 按照min_distance原则，选择最佳的Harris角点filtered_coords = []for i in index:if allowed_location[coords[i, 0], coords[i, 1]] == 1:filtered_coords.append(coords[i])allowed_location[(coords[i, 0] - min_dist):(coords[i, 0] + min_dist),(coords[i, 1] - min_dist):(coords[i, 1] + min_dist)] = 0return filtered_coords# 绘制图像中检测到的角点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()if "__main__" == __name__:# 打开灰度图convert('L')im = array(Image.open("jmuIma/zzw0.jpg").convert('L'))harrisim = compute_harris_response(im)filtered_coords = get_harris_points(harrisim, 6)plot_harris_points(im, filtered_coords)def get_descriptors(image, filtered_coords, wid=5):"""返回点周围2*wid+1个像素的值"""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 descdef 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_valuendx = argsort(-d)matchscores = ndx[:, 0]return matchscoresdef 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]] != n:matches_12[n] = -1return matches_12def appendimages(img1, img2):"""返回两幅图像并拼接成一幅新图像"""rows1 = img1.shape[0]rows2 = img2.shape[0]if rows1 < rows2:img1 = concatenate((img1, zeros((rows2 - rows1, img1.shape[1]))), axis=0)elif rows1 > rows2:img2 = concatenate((img2, zeros((rows1 - rows2, img2.shape[1]))), axis=0)return concatenate((img1, img2), axis=1)def plot_matches(img1, img2, locs1, locs2, matchscores, show_below=True):"""显示一幅带有连接匹配之间连线的图片"""img3 = appendimages(img1, img2)if show_below:img3 = vstack((img3, img3))imshow(img3)cols1 = img1.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')if "__main__" == __name__:# 打开灰度图convert('L')img1 = array(Image.open("../jmuIma/zzw0.jpg").convert('L'))img2 = array(Image.open("../jmuIma/zzw1.jpg").convert('L'))wid = 5harrisim = compute_harris_response(img1, 5)filtered_coords1 = get_harris_points(harrisim, wid + 1)d1 = get_descriptors(img1, filtered_coords1, wid)harrisim = compute_harris_response(img2, 5)filtered_coords2 = get_harris_points(harrisim, wid + 1)d2 = get_descriptors(img2, filtered_coords2, wid)print("starting matching")matches = match_twosided(d1, d2)figure()gray()plot_matches(img1, img2, filtered_coords1, filtered_coords2, matches)show()

效果