Harris角点检测、Sift图像特征匹配、地理特征匹配
Harris角点检测
什么是角点
下面有两幅不同视角的图像，通过找出对应的角点进行匹配。
角点在保留图像图形重要特征的同时，可以有效地减少信息的数据量，使其信息的含量很高，有效地提高了计算的速度，有利于图像的可靠匹配，使得实时处理成为可能。
我们可以直观的概括下角点所具有的特征：
- 轮廓之间的交点；
- 对于同一场景，即使视角发生变化，通常具备稳定性质的特征；
- 该点附近区域的像素点无论在梯度方向上还是其梯度幅值上有着较大变化；
角点检测算法基本思想
算法基本思想是使用一个固定窗口在图像上进行任意方向上的滑动，比较滑动前与滑动后两种情况，窗口中的像素灰度变化程度，如果存在任意方向上的滑动，都有着较大灰度变化，那么我们可以认为该窗口中存在角点。
角点检测算法基本原理
当窗口发生[u,v]移动时，那么滑动前与滑动后对应的窗口中的像素点灰度变化描述如下：
[u,v]是窗口的偏移量
(x,y)是窗口内所对应的像素坐标位置，窗口有多大，就有多少个位置
w(x,y)是窗口函数，最简单情形就是窗口内的所有像素所对应的w权重系数均为1。但有时候，我们会将w(x,y)函数设定为以窗口中心为原点的二元正态分布。如果窗口中心点是角点时，移动前与移动后，该点的灰度变化应该最为剧烈，所以该点权重系数可以设定大些，表示窗口移动时，该点在灰度变化贡献较大；而离窗口中心(角点)较远的点，这些点的灰度变化几*缓，这些点的权重系数，可以设定小点，以示该点对灰度变化贡献较小，那么我们自然想到使用二元高斯函数来表示窗口函数，所以通常窗口函数有如下两种形式：
根据上述表达式，当窗口处在平坦区域上滑动，可以想象的到，灰度不会发生变化，那么E(u,v) = 0；如果窗口处在比纹理比较丰富的区域上滑动，那么灰度变化会很大。算法最终思想就是计算灰度发生较大变化时所对应的位置，当然这个较大是指针任意方向上的滑动，并非单指某个方向。
I(x+u,y+v)泰勒展开可得：
当发生微小位移时，忽略无穷小量，写成矩阵形式：
所以E(u,v)表达式可以更新为：
矩阵M为:
M分解可得：
可以得出下列结论：
特征值都比较大时，即窗口中含有角点
特征值一个较大，一个较小，窗口中含有边缘
特征值都比较小，窗口处在平坦区域
可以通过判断R的值来判断某个点是不是角点了。
角点：R为大数值整数
边缘：R为大数值负数
平坦区：绝对值R是小数值
Harris角点检测代码
# -*- coding: utf-8 -*-
from pylab import *
from PIL import Image
from PCV.localdescriptors import harris
# 读入图像
im = array(Image.open('photo/photo1.jpg').convert('L'))
# 检测harris角点
harrisim = harris.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 = harris.get_harris_points(harrisim, 6, thres)
subplot(1, 4, i+2)
imshow(im)
print(im.shape)
plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '*')
axis('off')
show()
原图及Harris角点检测图
原图：
Harris角点检测图：
结果分析
实验过程中，选取像素值高于阈值的所有图像点作为角点，不过角点周围也有很多响应值很大的角点，使得角点比较密集，角点周围的点的角点响应函数值一般也会比较大，这样角点周围也会有很多点被当作是角点，使得检测的角点很密集。可以采用非极大值抑制的方法去除一些多余的角点，从而使更有特征的角点显示出来。
图像特征匹配
传统图像处理中图像特征匹配有三个基本步骤：特征提取、特征描述和特征匹配。特征提取就是从图像中提取出关键点（或特征点、角点）等。特征描述就是用一组数学向量对特征点进行描述，其主要保证不同的向量和不同的特征点之间是一种对应的关系，同时相似的关键点之间的差异尽可能小。特征匹配其实就是特征向量之间的距离计算，常用的距离有欧氏距离、汉明距离、余弦距离等。SIFT算法又叫尺度不变特征变换匹配算法， SIFT特征对于旋转和尺度均具有不变性，并且对于噪声、视角变化和光照变化具有良好的鲁棒性。
SIFT算法原理
①尺度空间极值检测——即找特征点
对于二维图像***I***(x，y)，建立图像的DOG （difference of guassians，即高斯差分）金字塔，DOG尺度空间含义为：可以用一个尺度空间的高斯函数和图像的卷积来表示。
G为尺度可变的高斯函数，***I***为空间坐标，其中西伽马为尺度。
为确定特征点所在的位置，首先需要建立一个高斯金字塔。具体可以参考文献。
得到高斯金字塔后，再通过两个相邻的高斯尺度空间做差，得到高斯差分DOC金字塔，整个过程的公式如下:
高斯差分金字塔建立后，特征点就是DOG尺度空间众多极值点，（个人认为对其进行求偏导数即可得到特征点的位置坐标），查找该极值点需要把每个点与邻域内的周围的26个点进行比较，这些点包括，在同一尺度上相邻的8个点，以及相邻尺度上相邻的18个点。
②特征点的尺度方向的确定
得到了特征点的坐标是完全不够的，必须要增加方向尺度信息。
- 采用有限差分的方法，求出在以特征点为圆心，以3倍西伽马为半径的范围内的图像梯度的幅值和相位。
- 利用直方图统计方法，求助邻域内所有像素点的梯度方向以及幅值。特征点的主方向就是直方图的峰值所代表的方向，确定了主方向就可以使SIFT算法具备旋转不变性。
- 主方向的计算公式为：（注：L表示特征点的尺度）
③特征向量的生成
- 特征向量最终是通过求得的特征点的邻域梯度信息来计算的。
- 先把坐标轴位置旋转到特征点所在的主方向上
- 接着以特征点为圆心，选择特征点附近的16个点作为种子点，分别求出8个方向上的梯度大小
- 最后得到的128维向量即为所求的特征向量
④特征点的匹配
通常采用最邻近的方法，即查找每一个特征点在另外一副图像中的最邻近。即最短的欧式距离。欧式距离的含义就是：两点连线的距离长度，具体如下：
- 对目标图像A中的某个特征点，求出该点与图像B中的所有特征点的特征向量的欧式距离；
- 将所求得的欧式距离值的大小进行排序；
- 找出最小的和次小(倒数第二小的)欧式距离值对于的目标图B中的特征点，并对这两个距离的比值进行计算；
- 假设该比值小于某个阈值，则这两个点为匹配点，否则不匹配。
（注：该阈值是一个经验阈值：是影响误匹配的主要原因，实验中一般选取的阈值大小是0.6左右）
SIFT图像匹配代码
from PIL import Image
from pylab import *
import sys
from PCV.localdescriptors import sift
if len(sys.argv) >= 3:
im1f, im2f = sys.argv[1], sys.argv[2]
else:
im1f = 'photo/photo8.jpg'
im2f = 'photo/photo4.jpg'
im1 = array(Image.open(im1f))
im2 = array(Image.open(im2f))
sift.process_image(im1f, 'out_sift_1.txt')
l1, d1 = sift.read_features_from_file('out_sift_1.txt')
figure()
gray()
subplot(121)
sift.plot_features(im1, l1, circle=False)
sift.process_image(im2f, 'out_sift_2.txt')
l2, d2 = sift.read_features_from_file('out_sift_2.txt')
subplot(122)
sift.plot_features(im2, l2, circle=False)
#matches = sift.match(d1, d2)
matches = sift.match_twosided(d1, d2)
print('{} matches'.format(len(matches.nonzero()[0])))
figure()
gray()
sift.plot_matches(im1, im2, l1, l2, matches, show_below=True)
show()
SIFT图像匹配测试结果
地理特征匹配
地理特征匹配代码
from numpy import zeros
from pylab import *
from PIL import Image
from PCV.localdescriptors import sift
from PCV.tools import imtools
import pydot
import os
os.environ["PATH"] += os.pathsep + 'D:/Graphviz/bin/'
""" This is the example graph illustration of matching images from Figure 2-10.
To download the images, see ch2_download_panoramio.py."""
download_path = "C:\\PycharmProjects\\untitled\\photo" # set this to the path where you downloaded the panoramio images
path = "C:\\PycharmProjects\\untitled\\photo\\" # path to save thumbnails (pydot needs the full system path)
# list of downloaded filenames
imlist = imtools.get_imlist(download_path)
nbr_images = len(imlist)
# extract features
featlist = [imname[:-3] + 'sift' for imname in imlist]
for i, imname in enumerate(imlist):
sift.process_image(imname, featlist[i])
matchscores = zeros((nbr_images, nbr_images))
for i in range(nbr_images):
for j in range(i, nbr_images): # only compute upper triangle
print('comparing ', imlist[i], imlist[j])
l1, d1 = sift.read_features_from_file(featlist[i])
l2, d2 = sift.read_features_from_file(featlist[j])
matches = sift.match_twosided(d1, d2)
nbr_matches = sum(matches > 0)
print('number of matches = ', nbr_matches)
matchscores[i, j] = nbr_matches
print("The match scores is: \n", matchscores)
# copy values
for i in range(nbr_images):
for j in range(i + 1, nbr_images): # no need to copy diagonal
matchscores[j, i] = matchscores[i, j]
#可视化
threshold = 2 # min number of matches needed to create link
g = pydot.Dot(graph_type='graph') # don't want the default directed graph
for i in range(nbr_images):
for j in range(i + 1, nbr_images):
if matchscores[i, j] > threshold:
# first image in pair
im = Image.open(imlist[i])
im.thumbnail((100, 100))
filename = path + str(i) + '.png'
im.save(filename) # need temporary files of the right size
g.add_node(pydot.Node(str(i), fontcolor='transparent', shape='rectangle', image=filename))
# second image in pair
im = Image.open(imlist[j])
im.thumbnail((100, 100))
filename = path + str(j) + '.png'
im.save(filename) # need temporary files of the right size
g.add_node(pydot.Node(str(j), fontcolor='transparent', shape='rectangle', image=filename))
g.add_edge(pydot.Edge(str(i), str(j)))
g.write_png('jmu.png')
测试图片及测试结果
测试原图：
测试结果：
结果分析
图左部分，建筑物在远处和近处能识别出是同一个建筑物，说明sift在一定程度上可以减少距离对图片匹配的影响。中间部分建筑从左侧和中间远近拍摄均能识别是同一个建筑，说明采用henssian矩阵获取图像局部最值还是十分稳定的。