SIFT特征提取实现地理位置的标注

1 SIFT特征提取分析

SIFT（Scale-invariant feature transform)是一种检测局部特征的算法，该算法通过求一幅图中的特征点（interest points,or corner points）及其有关scale 和 orientation 的描述子得到特征并进行图像特征点匹配，获得了良好效果。

算法描述：

SIFT特征不只具有尺度不变性，即使改变旋转角度，图像亮度或拍摄视角，仍然能够得到好的检测效果。整个算法分为以下几个部分：

1.构建尺度空间

这是一个初始化操作，尺度空间理论目的是模拟图像数据的多尺度特征。

高斯卷积核是实现尺度变换的唯一线性核，于是一副二维图像的尺度空间定义为：

$L(x,y,\sigma)=G(x,y,\sigma)\ast I(x, y)$

其中 $G(x,y,\sigma)$ 是尺度可变的高斯函数： $G(x,y,\sigma)=\frac{1}{2\pi\sigma^{2}}e^{-(x^{2}+y^{2})/2\sigma^{2}}$ 。

$(x, y)$ 是空间坐标。 $\sigma$ 大小决定图像的平滑程度，大尺度对应图像的盖帽特征，小尺度对应图像的细节特征。大的 $\sigma$ 值对应粗糙尺度（低分辨率），反之，对应精细尺度（高分辨率）。为了有效的在尺度空间检测到稳定的关键点，提出了高斯差分尺度空间（DOG scale-space)。利用不同尺度的高斯差分核与图像卷积生成。公式如下：

$D(x,y,\sigma)=(G(x,y,k\sigma)-G(x,y,\sigma)) \ast I(x,y) =L(x,y,k\sigma) - L(x,y,\sigma)$

关于尺度空间的理解说明：2kσ中的2是必须的，尺度空间是连续的。在 Lowe的论文中，将第0层的初始尺度定为1.6（最模糊），图片的初始尺度定为0.5（最清晰）. 在检测极值点前对原始图像的高斯平滑以致图像丢失高频信息，所以 Lowe 建议在建立尺度空间前首先对原始图像长宽扩展一倍，以保留原始图像信息，增加特征点数量。尺度越大图像越模糊。下图所示不同σ下图像尺度空间：

图像金字塔的建立：对于一幅图像I，建立其在不同尺度的图像，也称为子八度（octave），这是为了scale-invariant，也就是在任何尺度都能够有对应的特征点，第一个子八度的scale为原图大小，后面每个octave为上一个octave降采样的结果，即为原图的四分之一（长宽分别减半），构成下一个子八度（高一层金字塔）。

尺度空间的所有取值， $2^{i-1}(\sigma,k\sigma, k\sigma^{2},...k^{n-1}\sigma),k=2^{1/s}$ ,i为octave的塔数（第几个塔），s为每塔层数，由图片size决定建几个塔，每塔几层图像(S一般为3-5层)。0塔的第0层是原始图像(或你double后的图像)，往上每一层是对其下一层进行Laplacian变换（高斯卷积，其中σ值渐大，例如可以是σ, k*σ, k*k*σ…），直观上看来越往上图片越模糊。塔间的图片是降采样关系，例如1塔的第0层可以由0塔的第3层down sample得到，然后进行与0塔类似的高斯卷积操作。

2.LOG近似DOG找到关键点<检测DOG尺度空间极值点>

为了寻找尺度空间的极值点，每个采样点要和它所有的相邻点比较，看其是否比它的图像域和尺度域的相邻点大或者小。如图所示，中间的检测点和它同尺度的8个相邻点和上下相邻尺寸对应的9x2个点共同26个点比较，以确保在尺度空间和二维图像空间都检测到极值点。一个点如果在DOG尺度空间本层以及上下两层的26个领域中是最大或最小值时，就认为该点是图像在该尺度下的一个特征点,如图所示。

同一组中的相邻尺度（由于k的取值关系，肯定是上下层）之间进行寻找

在极值比较的过程中，每一组图像的首末两层是无法进行极值比较的，为了满足尺度变化的连续性。

3.除去不好的特征点

这一步本质上是要去掉DOG局部曲率非常不对称的像素。通过拟合三维二次函数以精确确定关键点的位置和尺寸（达到亚像素精度），同时去除低对比度的关键点和不稳定的边缘响应点（因为DOG算子会产生较强的边缘响应），以增强匹配稳定性、提高抗噪声能力，在这里我们使用近似Harris Corner检测器。

①空间尺度函数泰勒展开如下 $D(x) = D + \frac{\partial{D}^{T}}{\partial x}x + \frac{1}{2}x^{T}\frac{\partial ^{2}D}{\partial x^{2}}x$ ,对上式求导，并令其为0，得到精确的位置，得 $\hat{x}=-\frac{\partial ^{2}D^{-1}}{\partial x^{2}}\frac{\partial D}{\partial x}$ .

②在已经检测到的特征点中，要去掉低对比度的特征点和不稳定的边缘响应点。去除低对比度的点：把公式（2）代入公式（1）。即在DOG Space的极值点处D(x)取值，只取前两项可得：

$D(\hat x)=D+\frac{1}{2} \frac{\partial D^{T}}{\partial x} \hat x$ .

如果 $\mid D(\hat x)\mid \geq 0.03$ ,该特征点就保留下来，否则丢弃。

③边缘响应的去除

一个定义不好的高斯差分算子的极值在横跨边缘的地方有较大的主曲率，而在垂直边缘的方向有较小的主曲率。主曲率通过一个2x2的hessian矩阵H求出：

$H = \begin{bmatrix} D_{xx} & D_{xy}\\ D_{xy}&D_{yy} \end{bmatrix}$ ，导数由采样点相邻差估计得到。

D的主曲率和H的特征值成正比，令 $\alpha$ 为较大特征值， $\beta$ 为较小的特征值，则

$Tr(H) = D_{xx} + D_{yy} = \alpha + \beta$ ,

$Det(H)=D_{xx}D_{yy} - (D_{xy})^{2} = \alpha \beta$ .

令 $\alpha = \gamma \beta$ ,则：

$\frac{Tr(H)^{2}}{Det(H)} = \frac{(\alpha+\beta)^2}{\alpha \beta} = \frac{(r\beta+\beta)^{2}}{r\beta^{2}}=\frac{(r+1)^{2}}{r}$

$\frac{(r+1)^{2}}{r}$ 的值在两个特征值相等的时候最小，随着r的增大而增大，因此为了检测主曲率是否在某阈值r下，只需检测 $\frac{Tr(H)^{2}}{Det(H)} < \frac{(r+1)^{2}}{r}$ .

4. 给特征点赋值一个128维方向参数

上一步确定了每一幅图中的特征点，为每个特征点计算一个方向，依照这个方向做进一步的计算，利用关键点领域像素的梯度方向分布特征为每个关键点指定方向参数，使算子具备旋转不变性。

$m(x,y)=\sqrt{(L(x+1, y)-L(x-1, y))^{2}+(L(x,y+1)-L(x,y-1))^{2}}$

$\theta(x,y)=\alpha \tan2((L(x,y+1)-L(x,y-1))/(L(x+1, y)-L(x-1, y)))$

为(x,y)处梯度的模值和方向公式。其中 $\angle$ 所用的尺度为每个关键点各自所在的尺度。至此，图像的关键点已经检测完毕，每个关键点有三个信息：位置、所处尺度、方向，由此可以确定一个SIFT特征区域。

梯度直方图的范围是0～360度，其中每10度一个柱，总共36个柱。随着距中心点越远的领域其对直方图的贡献也响应减小.Lowe论文中还提到要使用高斯函数对直方图进行平滑，减少突变的影响。

实际计算时，我们在以关键点为中心的邻域窗口内采样，并用直方图统计邻域像素的梯度方向。梯度直方图的范围是0～360度，其中每45度一个柱，总共8个柱, 或者每10度一个柱，总共36个柱。Lowe论文中还提到要使用高斯函数对直方图进行平滑，减少突变的影响。直方图的峰值则代表了该关键点处邻域梯度的主方向，即作为该关键点的方向。

直方图中的峰值就是主方向，其他的达到最大值80%的方向可作为辅助方向。

由梯度方向直方图确定主梯度方向，该步中将建立所有scale中特征点的描述子(128维）。

通过对关键点周围图像区域分块，计算块内梯度直方图，生成具有独特性的向量，这个向量是该区域图像信息的一种抽象，具有唯一性。

5.关键点描述子的生成

首先将坐标轴旋转为关键点的方向，以确保旋转不变性。以关键点为中心取8x8的窗口。

16*16的图中其中1/4的特征点梯度方向及scale，右图为其加权到8个主方向后的效果。图左部分的中央为当前关键点的位置，每个小格代表关键点领域所在尺度空间的一个像素，利用公式求得每个像素的梯度幅值与梯度方向，箭头方向代表该像素的梯度方向，箭头长度代表梯度的模值，然后用高斯窗口对其进行加权运算。

图中蓝色的圈代表高斯加权的范围（越靠近关键点的像素梯度方向信息贡献越大）。然后在每4×4的小块上计算8个方向的梯度方向直方图，绘制每个梯度方向的累加值，即可形成一个种子点，如图右部分示。此图中一个关键点由2×2共4个种子点组成，每个种子点有8个方向向量信息。这种邻域方向性信息联合的思想增强了算法抗噪声的能力，同时对于含有定位误差的特征匹配也提供了较好的容错性。

计算keypoint周围的16*16的window中每一个像素的梯度，而且使用高斯下降函数降低远离中心的权重。

在每个4*4的1/16象限中，通过加权梯度值加到直方图8个方向区间中的一个，计算出一个梯度方向直方图。这样就可以对每个feature形成一个4*4*8=128维的描述子，每一维都可以表示4*4个格子中一个的scale/orientation. 将这个向量归一化之后，就进一步去除了光照的影响。

6.根据SIFT进行Match

生成了A、B两幅图的描述子，(分别使K1*128维和K2*128维），就将两图中各个scale的描述子进行匹配，匹配上128维即可表示两个特征点match上了。

当两幅图像的SIFT特征向量生成后，下一步我们采用关键点特征向量的欧式距离来作为两幅图像中关键点的相似性判定度量。取图像1中的某个关键点，并找出其与图像2中欧式距离最近的前两个关键点，在这两个关键点中，如果最近的距离除以次近的距离少于某个比例阈值，则接受这一对匹配点。降低这个比例阈值，SIFT匹配点数目会减少，但更加稳定。

2 地理位置坐标匹配

我们首先准备好15张建筑物数据集，然后用前面用到的SIFT描述子。

采用pydot工具包对SIFT匹配结果进行可视化：

from pylab import *
from PIL import Image
from PCV.localdescriptors import sift
from PCV.tools import imtools
import pydotdownload_path = "img3"  # set this to the path where you downloaded the panoramio images
path = "img3/"  # 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 triangleprint('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 diagonalmatchscores[j, i] = matchscores[i, j]# 可视化threshold = 2  # min number of matches needed to create linkg = pydot.Dot(graph_type='graph')  # don't want the default directed graphfor i in range(nbr_images):for j in range(i + 1, nbr_images):if matchscores[i, j] > threshold:# first image in pairim = Image.open(imlist[i])im.thumbnail((100, 100))    # 图片按比例缩放filename = path + str(i) + '.png'im.save(filename)  # need temporary files of the right sizeg.add_node(pydot.Node(str(i), fontcolor='transparent', shape='rectangle', image=filename))# second image in pairim = Image.open(imlist[j])im.thumbnail((100, 100))filename = path + str(j) + '.png'im.save(filename)  # need temporary files of the right size# im.show()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('JiMei_University.png')

效果图: