大家好，又见面了，我是你们的朋友全栈君。

传统图像处理中图像特征匹配有三个基本步骤：特征提取、特征描述和特征匹配。特征提取就是从图像中提取出关键点（或特征点、角点）等。特征描述就是用一组数学向量对特征点进行描述，其主要保证不同的向量和不同的特征点之间是一种对应的关系，同时相似的关键点之间的差异尽可能小。特征匹配其实就是特征向量之间的距离计算，常用的距离有欧氏距离、汉明距离、余弦距离等。

SIFT算法又叫尺度不变特征变换匹配算法， SIFT特征对于旋转和尺度均具有不变性，并且对于噪声、视角变化和光照变化具有良好的鲁棒性，所以我们今天来学习一下SIFT算法。

二、SIFT算法原理

四、实现

1.vlfeat安装

3.关于图像

一、SIFT简介

SIFT（Scale Invariant Feature Transform，尺度不变特征变换匹配算法）是由David G. Lowe教授在1999年（《Object Recognition from Local Scale-Invariant Features》）提出的高效区域检测算法，在2004年（《Distinctive Image Features from Scale-Invariant Keypoints》）得以完善。

SIFT可以应用到物体辨识、机器人地图感知与导航、影像缝合、3D模型建立、手势辨识、影像追踪和动作比对等方向。

SIFT算法的特点:

稳定性
独特性
多量性
高速性
可扩展性

SIFT算法可以的解决问题：

目标的旋转、缩放、平移（RST）
图像放射/投影变换（视点viewpoint）
光照影响（illumination）
部分目标遮挡（occlusion）
杂物场景（clutter）
噪声

二、SIFT算法原理

1.检测尺度空间极值

检测尺度空间极值就是搜索所有尺度上的图像位置，通过高斯微分函数来识别对于尺度和旋转不变的兴趣点。其主要步骤可以分为建立高斯金字塔、生成DOG高斯差分金字塔和DOG局部极值点检测。为了让大家更清楚，我先简单介绍一下尺度空间，再介绍主要步骤。

（1）尺度空间

一个图像的尺度空间，定义为一个变化尺度的高斯函数与原图像的卷积。即：

，

其中，*表示卷积计算。

其中，m、n表示高斯模版的维度，(x,y)代表图像像素的位置。为尺度空间因子，值越小表示图像被平滑的越少，相应的尺度就越小。小尺度对应于图像的细节特征，大尺度对应于图像的概貌特征，效果如下图所示，尺度从左到右，从上到下，一次增大。

（2）建立高斯金字塔

尺度空间在实现时，使用高斯金字塔表示，高斯金字塔的构建分为两部分：

1.对图像做不同尺度的高斯模糊

2.对图像做降采样（隔点采样）

图像的金字塔模型是指，将原始图像不断降阶采样，得到一系列大小不一的图像，由大到小，从下到上构成的塔状模型。原图像为金子塔的第一层，每次降采样所得到的新图像为金字塔的上一层(每层一张图像)，每个金字塔共n层。金字塔的层数根据图像的原始大小和塔顶图像的大小共同决定。

为了让尺度体现其连续性，高斯金字塔在简单降采样的基础上加上了高斯滤波。如上图所示，将图像金字塔每层的一张图像使用不同参数做高斯模糊，使得金字塔的每层含有多张高斯模糊图像，将金字塔每层多张图像合称为一组(Octave)，金字塔每层只有一组图像，组数和金字塔层数相等，每组含有多层Interval图像。

高斯图像金字塔共o组、s层，则有：

其中，σ表示尺度空间坐标，s表示sub-level层坐标，表示初始尺度，S表示每组层数（一般为3~5）

（3）建立DOG高斯差分金字塔

为了有效提取稳定的关键点，利用不同尺度的高斯差分核与卷积生成。

DOG函数：

DOG在计算上只需相邻高斯平滑后图像相减，因此简化了计算！

可以通过高斯差分图像看出图像上的像素值变化情况。（如果没有变化，也就没有特征。特征必须是变化尽可能多的点。）DOG图像描绘的是目标的轮廓。

（4）DOG局部极值检测

特征点是由DOG空间的局部极值点组成的。为了寻找DOG函数的极值点，每一个像素点要和它所有的相邻点比较，看其是否比它的图像域和尺度域的相邻点大或者小。

中间的检测点和它同尺度的8个相邻点和上下相邻尺度对应的9×2个点共26个点比较，以确保在尺度空间和二维图像空间都检测到极值点。

2.关键点的精确定位

以上方法检测到的极值点是离散空间的极值点，以下通过拟合三维二次函数来精确确定关键点的位置和尺度，同时去除低对比度的关键点和不稳定的边缘响应点(因为DOG算子会产生较强的边缘响应)，以增强匹配稳定性、提高抗噪声能力。

（1）关键点的精确定位

利用已知的离散空间点插值得到的连续空间极值点的方法叫做子像素插值（Sub-pixel Interpolation）。

为了提高关键点的稳定性，需要对尺度空间DOG函数进行曲线拟合。利用DOG函数在尺度空间的Taylor展开式(拟合函数)为：

其中，。求导并让方程等于零，可以得到极值点的偏移量为：

对应极值点，方程的值为：

其中, 代表相对插值中心的偏移量，当它在任一维度上的偏移量大于0.5时（即x或y或），意味着插值中心已经偏移到它的邻近点上，所以必须改变当前关键点的位置。同时在新的位置上反复插值直到收敛；也有可能超出所设定的迭代次数或者超出图像边界的范围，此时这样的点应该删除，在Lowe中进行了5次迭代。另外，过小的点易受噪声的干扰而变得不稳定，所以将小于某个经验值(Lowe论文中使用0.03，Rob Hess等人实现时使用0.04/S)的极值点删除。同时，在此过程中获取特征点的精确位置(原位置加上拟合的偏移量)以及尺度()。

（2）去除边缘响应

由于DOG函数在图像边缘有较强的边缘响应，因此需要排除边缘响应DOG函数的峰值点在边缘方向有较大的主曲率，而在垂直边缘的方向有较小的主曲率。主曲率可以通过计算在该点位置尺度的2×2的Hessian矩阵得到，导数由采样点相邻差来估计：

表示DOG金字塔中某一尺度的图像x方向求导两次。

D的主曲率和H的特征值成正比。令 α ，β为特征值，则

该值在两特征值相等时达最小。Lowe论文中建议阈值T为1.2，即时保留关键点，反之剔除。

在Lowe的论文中，取r＝10。下图右侧为消除边缘响应后的关键点分布图。

3.关键点主方向分配

关键点主方向分配就是基于图像局部的梯度方向，分配给每个关键点位置一个或多个方向。所有后面的对图像数据的操作都相对于关键点的方向、尺度和位置进行变换，使得描述符具有旋转不变性。

对于在DOG金字塔中检测出的关键点，采集其所在高斯金字塔图像3σ邻域窗口内像素的梯度和方向分布特征。梯度的模值和方向如下：

L为关键点所在的尺度空间值，按Lowe的建议，梯度的模值m(x,y)按的高斯分布加成，按尺度采样的3σ原则，邻域窗口半径为。

在完成关键点的梯度计算后，使用直方图统计邻域内像素的梯度和方向。梯度直方图将0~360度的方向范围分为36个柱(bins)，其中每柱10度。如下图所示，直方图的峰值方向代表了关键点的主方向，(为简化，图中只画了八个方向的直方图)。

方向直方图的峰值代表了该特征点处邻域梯度的方向，以直方图中最大值作为该关键点的主方向。为了增强匹配的鲁棒性，只保留峰值大于主方向峰值80％的方向作为该关键点的辅方向。Lowe的论文指出大概有15%关键点具有多方向，但这些点对匹配的稳定性至为关键。检测结果如下图：

至此，将检测出的含有位置、尺度和方向的关键点即是该图像的SIFT特征点。

4.关键点的特征描述

通过以上步骤，对于每一个关键点，拥有三个信息：位置、尺度以及方向。接下来就是为每个关键点建立一个描述符，用一组向量将这个关键点描述出来，使其不随各种变化而改变，比如光照变化、视角变化等。这个描述子不但包括关键点，也包含关键点周围对其有贡献的像素点，并且描述符应该有较高的独特性，以便于提高特征点正确匹配的概率。

SIFT描述子是关键点邻域高斯图像梯度统计结果的一种表示。通过对关键点周围图像区域分块，计算块内梯度直方图，生成具有独特性的向量，这个向量是该区域图像信息的一种抽象，具有唯一性。

Lowe建议描述子使用在关键点尺度空间内4*4的窗口中计算的8个方向的梯度信息，共4*4*8=128维向量表征。表示步骤如下：

（1）计算描述子所需的图像区域

特征描述子与特征点所在的尺度有关，因此，对梯度的求取应在特征点对应的高斯图像上进行。将关键点附近的邻域划分为d*d(Lowe建议d=4)个子区域，每个子区域做为一个种子点，每个种子点有8个方向。每个子区域的大小与关键点方向分配时相同，即每个区域有个子像素，为每个子区域分配边长为的矩形区域进行采样(个子像素实际用边长为矩形区域进行采样，考虑到实际计算时，需要采用双线性插值，所需图像窗口边长为。在考虑到旋转因素(方便下一步将坐标轴旋转到关键点的方向)，如下图6.1所示，实际计算所需的图像区域半径为：

计算结果四舍五入取整。

（2）将坐标轴旋转为关键点的方向

将坐标轴旋转为关键点的方向是为了确保旋转不变性，如下图所示：

旋转后邻域内采样点的新坐标为：

（3）将邻域内的采样点分配到对应的子区域内

将子区域内的梯度值分配到8个方向上，计算其权值。旋转后的采样点坐标在半径为radius的圆内被分配到的子区域，计算影响子区域的采样点的梯度和方向，分配到8个方向上。

旋转后的采样点落在子区域的下标为

(6-3)

Lowe建议子区域的像素的梯度大小按的高斯加权计算，即

(6-4)

其中a，b为关键点在高斯金字塔图像中的位置坐标。

（4）插值计算每个种子八个方向的梯度。

如上图所示，将所得采样点在子区域中的下标(x”,y”)(图中蓝色窗口内红色点)线性插值，计算其对每个种子点的贡献。如图中的红色点，落在第0 行和第1 行之间，对这两行都有贡献。对第0 行第3 列种子点的贡献因子为dr，对第1 行第3 列的贡献因子为1-dr，同理，对邻近两列的贡献因子为dc 和1-dc，对邻近两个方向的贡献因子为do 和1-do。则最终累加在每个方向上的梯度大小为：

其中，k,m,n为0或为1.

（5）描述符向量元素门限化

即把方向直方图每个方向上梯度幅值限制在一定门限值一下（门限一般取0.2）

（6）描述符向量元素归一化

特征向量形成后，为了去除光照变化的影响，需要对它们进行归一化处理，对于图像灰度值整体漂移，图像各点的梯度是邻域像素相减得到，所以也能去除。

得到的描述子向量：

归一化后的特征向量：

则

三、关键点匹配

分别对模板图（参考图，reference image）和实时图（观测图， observation image）建立关键点描述子集合。目标的识别是通过两点集内关键点描述子的比对来完成。具有128维的关键点描述子的相似性度量采用欧式距离。

模板图中关键点描述子：

实时图中关键点描述子：

任意两描述子相似性度量：

要得到配对的关键点描述子，需满足

关键点的匹配可以采用穷举法来完成，但是这样耗费的时间太多，一般都采用kd树的数据结构来完成搜索。搜索的内容是以目标图像的关键点为基准，搜索与目标图像的特征点最邻近的原图像特征点和次邻近的原图像特征点。

四、实现

1.检测感兴趣点

为了计算图像的SIFT特征，我们用开源工具包VLFeat。用Python重新实现SIFT特征提取的全过程不会很高效，而且也超出了本书的范围。VLFeat可以在www.vlfeat.org上下载，它的二进制文件可以用于一些主要的平台。这个库是用C写的，不过我们可以利用它的命令行接口。此外，它还有Matlab接口。

# -*- coding: utf-8 -*-
from PIL import Image
from pylab import *
from PCV.localdescriptors import sift
from PCV.localdescriptors import harris

# 添加中文字体支持
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)

imname = 'D:/test/111.jpg'
im = array(Image.open(imname).convert('L'))
sift.process_image(imname, '111.sift')
l1, d1 = sift.read_features_from_file('111.sift')

figure()
gray()
subplot(131)
sift.plot_features(im, l1, circle=False)
title(u'SIFT特征',fontproperties=font)
subplot(132)
sift.plot_features(im, l1, circle=True)
title(u'用圆圈表示SIFT特征尺度',fontproperties=font)

# 检测harris角点
harrisim = harris.compute_harris_response(im)

subplot(133)
filtered_coords = harris.get_harris_points(harrisim, 6, 0.1)
imshow(im)
plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '*')
axis('off')
title(u'Harris角点',fontproperties=font)

show()

为了将sift和Harris角点进行比较，将Harris角点检测的显示在了图像的最后侧。正如你所看到的，这两种算法选择了不同的坐标。由图可以看出，与Harris角点检测相比，sift提取出来的特征点信息更多，而且更加精准，这是因为sift的特征点提取步骤比Harris的步骤复杂的多，它需要建立高斯图像金字塔和高斯差分金字塔之后再检测极值，而Harris角点只是对原图进行角点检测和变化。

2.描述子匹配

from PIL import Image
from pylab import *
import sys
from PCV.localdescriptors import sift
from PIL import Image
from pylab import *
from numpy import *
import os

if len(sys.argv) >= 3:
  im1f, im2f = sys.argv[1], sys.argv[2]
else:
  im1f = 'D:/test/111.jpg'
  im2f = 'D:/test/222.jpg'

  #im1f = 'D:/test/change/21.jpg'
  #im2f = 'D:/test/change/22.jpg'
#  im1f = '../data/crans_1_small.jpg'
#  im2f = '../data/crans_2_small.jpg'
#  im1f = '../data/climbing_1_small.jpg'
#  im2f = '../data/climbing_2_small.jpg'
im1 = array(Image.open(im1f))
im2 = array(Image.open(im2f))

sift.process_image(im1f, 'out_sift_3.txt')
l1, d1 = sift.read_features_from_file('out_sift_3.txt')
figure()
gray()
subplot(121)
sift.plot_features(im1, l1, circle=False)

sift.process_image(im2f, 'out_sift_4.txt')
l2, d2 = sift.read_features_from_file('out_sift_4.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检测两张图片的感兴趣点的结果图：

由图可以看出，SIFT可以检测出较多的特征点。

下图左侧为SIFT算法匹配结果，右侧为Harris角点匹配结果：

由上面两张图对比可知，SIFT算法匹配出的特征点更多，这是因为SIFT算法具有尺度和旋转不变性，即使两张图大小不一样、角度不一致也不会影响匹配结果，而Harris角点对尺度变化非常敏感，当遇到尺度变化较大时，很多正确特征点无法检测出来。

3.地理标记图像匹配

我们可以对上面匹配后的图像进行连接可视化，要做到这样，我们需要在一个图中用边线表示它们之间是相连的。我们采用pydot工具包,它提供了GraphViz graphing库的Python接口。不要担心，它们安装起来很容易。

# -*- coding: utf-8 -*-
from pylab import *
from PIL import Image
from PCV.localdescriptors import sift
from PCV.tools import imtools
import pydot

""" This is the example graph illustration of matching images from Figure 2-10.
To download the images, see ch2_download_panoramio.py."""

#download_path = "panoimages"  # set this to the path where you downloaded the panoramio images
#path = "/FULLPATH/panoimages/"  # path to save thumbnails (pydot needs the full system path)

#download_path = "F:\\dropbox\\Dropbox\\translation\\pcv-notebook\\data\\panoimages"  # set this to the path where you downloaded the panoramio images
#path = "F:\\dropbox\\Dropbox\\translation\\pcv-notebook\\data\\panoimages\\"  # path to save thumbnails (pydot needs the full system path)
download_path = "D:/test/change"
path = "D:/test/change/"
# 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: %d", matchscores

#np.savetxt(("../data/panoimages/panoramio_matches.txt",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('change3.png')