python计算机视觉入门二：SIFT特征匹配及其实现
SIFT（尺度不变特征变换，Scale-Invariant Feature Transform）是在计算机视觉领域中检测和描述图像中局部特征的算法，该算法于1999年被David Lowe提出，并于2004年进行了补充和完善。该算法应用很广，如目标识别，自动导航，图像拼接，三维建模，手势识别，视频跟踪等。
一，sift算法的特征
1.SIFT特征是图像的局部特征，其对旋转、尺度缩放、亮度变化保持不变性，对视角变化、仿射变换、噪声也保持一定程度的稳定性；
2. 区分性（Distinctiveness）好，信息量丰富，适用于在海量特征数据库中进行快速、准确的匹配；
3. 多量性，即使少数的几个物体也可以产生大量的SIFT特征向量；
4.高速性，经优化的SIFT匹配算法甚至可以达到实时的要求；
5.可扩展性，可以很方便的与其他形式的特征向量进行联合。
二,SIFT算法的过程
1，提取关键点
尺度不变性------高斯金字塔
尺度空间主要思想是通过对原始图像进行尺度变换，获得图像多尺度下的空间表示。从而实现边缘，角点检测和不同分辨率上的特征提取，以满足特征点的尺度不变性。尺度空间中各尺度图像模糊程度逐渐变大，能够模拟人在距离目标由近到原时目标在视网膜上的形成过程。
尺度空间在实现高斯金字塔表示时，高斯金字塔的实现分为两个步骤：
1，对图像做不同尺度的高斯模糊
2，对图像做降采用
所谓图像金字塔化的构建：就是先进行图像平滑，再进行降采样，根据降采样率不同，所得到一系列尺寸逐渐减小的图像。
两者的区别：通过高斯模糊获取相同分辨率不同模糊程度的图像，从获取分辨率变化最大的点获得图像的轮廓特征，通过降采样获取不同尺寸的图像从而得到一系列图像尺寸变小的图像。
LOG（Laplassian of Gaussian）
通过上面的步骤建立的金字塔是具有不同尺寸，每一层只有一张图像的金字塔；通过结合两者有时得到LOG图像即高斯拉普拉斯变换图像。其步骤为：先将照片降采样，得到了不同分辨率下的图像金字塔。再对每层图像进行高斯卷积。这样一来，原本的图像金字塔每层只有一张图像，而卷积后，每层又增加了多张不同模糊程度下的照片。
差分金字塔
构造高斯差分图像的步骤是：在获得LOG图像后，用其相邻的图像进行相减，得到所有图像重新构造的金字塔就是DOG金字塔。
通过构建差分金字塔更好的获取关键点，进一步优化了LOG图像，减少计算量。
DOG局部极值点
当得到DOG金字塔后，我们接下来要做的是寻找DOG极值点。每个像素点与其周围的像素点比较，当其大于或者小于所有相邻点时，即为极值点。
极值点定位
通过DOG金字塔找到的极值点是属于离散空间上的点，不一定是真正意义上的极值点。我们可以利用已知的离散空间点插值得到的连续空间极值点的方法叫做子像素插值（Sub-pixel Interpolation），对尺度空间DoG函数进行曲线拟合，从而提高关键点的定位。
到这一步，得到的极值点是比较精确了，但不够准确。有些极值点不是我们想要的，当中就有一大部分是边缘区域产生的极值点。因为物体的边缘轮廓在灰度图中，存在着灰度值的突变例如图像边缘区域，这样的点在计算中就被“误以为”是特征值。因此，我们可以用Hessian矩阵来判断；
DOG函数的峰值点在边缘方向有较大的主曲率，而在垂直边缘的方向有较小的主曲率。主曲率可以通过计算在该点位置尺度的2X2的Hessian矩阵得到，导数由采样点相邻差来估计：
二，特征点描述
经过上述过程之后，基本上得到了我们想要的精确特征点了。接下来我们就要求它们的方向。
通过尺度不变性求极值点，可以使其具有缩放不变的性质。而利 用关键点邻域像素的梯度方向分布特性，可以为每个关键点指定方向参数 方向，从而使描述子对图像旋转具有不变性。
通过求每个极值点的梯度来为极值点赋予方向。
像素点的梯度表示
确定关键点的方向采用梯度直方图统计法，统计以关键点为原 点，一定区域内的图像像素点对关键点方向生成所作的贡献
关键点主方向：极值点周围区域梯度直方图的主峰值也是特征点方向 • 关键点辅方向：在梯度方向直方图中，当存在另一个相当于主峰值 80%能量的峰值时，则将这个方向认为是该关键点的辅方向。
3，关键点匹配
通过以上步骤，对于每一个关键点，拥有三个信息：位置、尺度以及方向。接下来就是为每个关键点建立一个描述符，使其不随各种变化而改变，比如光照变化、视角变化等等。并且描述符应该有较高的独特性，以便于提高特征点正确匹配的概率。 SIFI 描述子h(x, y,θ)是对特征点附近邻域内高斯图像梯度统计结果的一种表示，它是一个三维的阵列，但通常将它表示成一个矢量。矢量是通过对三维阵列按一定规律进行排列得到的。特征描述子与特征点所在的尺度有关，因此，对梯度的求取应在特征点对应的高斯图像上进行。
接下来就是生成特征匹配点。
在每子区域内计算8个方向的梯度方向直方图，绘制每个梯度方向的累加值，形成一个种子点。与求特征点主方向时有所不同，此时，每个子区域的梯度方向直方图将0°~360°划分为8个方向范围，每个范围为45°，这样，每个种子点共有8个方向的梯度强度信息。由于存在4X4（Bp X Bp）个子区域，所以，共有4X4X8=128个数据，最终形成128维的SIFT特征矢量。同样，对于特征矢量需要进行高斯加权处理，加权采用方差为mσBp/2的标准高斯函数，其中距离为各点相对于特征点的距离。使用高斯权重的是为了防止位置微小的变化给特征向量带来很大的改变，并且给远离特征点的点赋予较小的权重，以防止错误的匹配。
三，SIFT算法python代码
1，SIFT特征提取与Harris角点检测对比
from PIL import Image
from pylab import *
from numpy import *
import os
def process_image(imagename, resultname, params="--edge-thresh 10 --peak-thresh 5"):
""" 处理一幅图像，然后将结果保存在文件中"""
if imagename[-3:] != 'pgm':
#创建一个pgm文件
im = Image.open(imagename).convert('L')
im.save('tmp.pgm')
imagename ='tmp.pgm'
cmmd = str("D:\python计算机视觉\sift.exe " + imagename + " --output=" + resultname +
" " + params)
os.system(cmmd)
print('processed', imagename, 'to', resultname)
def read_features_from_file(filename):
"""读取特征属性值，然后将其以矩阵的形式返回"""
f = loadtxt(filename)
return f[:,:4], f[:,4:] #特征位置，描述子
def write_featrues_to_file(filename, locs, desc):
"""将特征位置和描述子保存到文件中"""
savetxt(filename, hstack((locs,desc)))
def plot_features(im, locs, circle=False):
"""显示带有特征的图像
输入：im（数组图像），locs（每个特征的行、列、尺度和朝向）"""
def draw_circle(c,r):
t = arange(0,1.01,.01)*2*pi
x = r*cos(t) + c[0]
y = r*sin(t) + c[1]
plot(x, y, 'b', linewidth=2)
imshow(im)
if circle:
for p in locs:
draw_circle(p[:2], p[2])
else:
plot(locs[:,0], locs[:,1], 'ob')
axis('off')
imname = 'empire.jpg'
im1 = array(Image.open(imname).convert('L'))
process_image(imname, 'empire.sift')
l1,d1 = read_features_from_file('empire.sift')
figure()
gray()
plot_features(im1, l1, circle=True)
show()
第一张图为用过Harris焦点检测提取的特征点，第二张图为通过SIFT特征提取获得的特征点 。
2，两张图像的SIFT特征匹配
sift.py
from PIL import Image
from pylab import *
from numpy import *
import os
def process_image(imagename, resultname, params="--edge-thresh 10 --peak-thresh 5"):
""" 处理一幅图像，然后将结果保存在文件中"""
if imagename[-3:] != 'pgm':
#创建一个pgm文件
im = Image.open(imagename).convert('L')
im.save('tmp.pgm')
imagename ='tmp.pgm'
cmmd = str("sift "+imagename+" --output="+resultname+" "+params)
os.system(cmmd)
print ('processed', imagename, 'to', resultname)
def read_features_from_file(filename):
"""读取特征属性值，然后将其以矩阵的形式返回"""
f = loadtxt(filename)
return f[:,:4], f[:,4:] #特征位置，描述子
def write_featrues_to_file(filename, locs, desc):
"""将特征位置和描述子保存到文件中"""
savetxt(filename, hstack((locs,desc)))
def plot_features(im, locs, circle=False):
"""显示带有特征的图像
输入：im（数组图像），locs（每个特征的行、列、尺度和朝向）"""
def draw_circle(c,r):
t = arange(0,1.01,.01)*2*pi
x = r*cos(t) + c[0]
y = r*sin(t) + c[1]
plot(x, y, 'b', linewidth=2)
imshow(im)
if circle:
for p in locs:
draw_circle(p[:2], p[2])
else:
plot(locs[:,0], locs[:,1], 'ob')
axis('off')
def match(desc1, desc2):
"""对于第一幅图像中的每个描述子，选取其在第二幅图像中的匹配
输入：desc1(第一幅图像中的描述子)，desc2(第二幅图像中的描述子)"""
desc1 = array([d/linalg.norm(d) for d in desc1])
desc2 = array([d/linalg.norm(d) for d in desc2])
dist_ratio = 0.6
desc1_size = desc1.shape
matchscores = zeros((desc1_size[0],1),'int')
desc2t = desc2.T #预先计算矩阵转置
for i in range(desc1_size[0]):
dotprods = dot(desc1[i,:],desc2t) #向量点乘
dotprods = 0.9999*dotprods
# 反余弦和反排序，返回第二幅图像中特征的索引
indx = argsort(arccos(dotprods))
#检查最近邻的角度是否小于dist_ratio乘以第二近邻的角度
if arccos(dotprods)[indx[0]] < dist_ratio * arccos(dotprods)[indx[1]]:
matchscores[i] = int(indx[0])
return matchscores
def match_twosided(desc1, desc2):
"""双向对称版本的match()"""
matches_12 = match(desc1, desc2)
matches_21 = match(desc2, desc1)
ndx_12 = matches_12.nonzero()[0]
# 去除不对称的匹配
for n in ndx_12:
if matches_21[int(matches_12[n])] != n:
matches_12[n] = 0
return matches_12
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 in range(len(matchscores)):
if matchscores[i]>0:
plot([locs1[i,0],locs2[matchscores[i,0],0]+cols1], [locs1[i,1],locs2[matchscores[i,0],1]],'c')
axis('off')
# 示例
im1f = 'crans_1_small.jpg'
im2f = 'crans_2_small.jpg'
im1 = array(Image.open(im1f))
im2 = array(Image.open(im2f))
process_image(im1f, 'out_sift_1.txt')
l1,d1 = read_features_from_file('out_sift_1.txt')
figure()
gray()
subplot(121)
plot_features(im1, l1, circle=False)
process_image(im2f, 'out_sift_2.txt')
l2,d2 = read_features_from_file('out_sift_2.txt')
subplot(122)
plot_features(im2, l2, circle=False)
matches = match_twosided(d1, d2)
print('{} matches'.format(len(matches.nonzero()[0])))
figure()
gray()
plot_matches(im1,im2,l1,l2,matches, show_below=True)
show()
3，SIFT特征做地理标记
# -*- coding: utf-8 -*-
from PIL import Image
from pylab import *
from numpy import *
import os
import pydotplus
import sift
def get_imlist(path):
return [os.path.join(path,f) for f in os.listdir(path) if f.endswith('.jpg')]
download_path = "D:/img"
path = "D:/img/"
# list of downloaded filenames
imlist = get_imlist(download_path)
nbr_images = len(imlist)
print(nbr_images)
# 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 craete link
g = pydotplus.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:
#图像对中的第一幅图像
im = Image.open(imlist[i])
im.thumbnail((100,100))
filename = path+str(i) + '.png'
im.save(filename) #需要一定大小的临时文件
g.add_node(pydotplus.Node(str(i), fontcolor='transparent',
shape='rectangle', image=filename))
#图像对中的第二幅图像
im = Image.open(imlist[j])
im.thumbnail((100,100))
filename = path+str(j) + '.png'
print(filename);
im.save(filename) #需要一定大小的临时文件
g.add_node(pydotplus.Node(str(j), fontcolor='transparent',
shape='rectangle', image=filename))
g.add_edge(pydotplus.Edge(str(i), str(j)))
g.write_png('result.png')
四，SIFT代码实现遇到的问题
1，Python实现SIFT匹配需要下载VFeat工具包，建议下载vlfeat-0.9.20-bin.tar.gz ，解压缩后，将vlfeat-0.9.20/bin/win64文件夹下的sift.exe和vl.dll拷贝到当前工作目录下。 我原来使用vlfeat-0.9.28-bin.tar.g出现生成的empire.sift文件为空的
2，在实现可视化连接图时，《python计算机视觉》下载pydot和GraphViz，但时由于我下载的pydot版本为1.41,运行时出现
find_GraphViz() is not defined 问题，原因是新版的pydot里已经不存在该函数了。解决方案：pip install pydotplust ;安装的时候会问要不要一起安装Graphviz.
3,在实现连接图出现，可以得出图片如下图，但是无法显示图片，控制台显示如下错误：
Warining: no or imprope .....原因是我用的是路径含有中文名，而且要注意：pydotplus的Node节点添加图片时，图片的路径需要为绝对路径，且分隔符为/。
参考资料：
https://blog.csdn.net/zhuxiaoyang2000/article/details/53930610
https://blog.csdn.net/weixin_38404120/article/details/73740612
《python计算机视觉》