Python SIFT特征匹配算法、匹配地理标记图像

一、SIFT算法的原理

1、检测适度空间极值

　　检测尺度空间极值就是搜索所有尺度上的图像位置，通过高斯微分函数来识别对于尺寸和旋转不变的兴趣点。其主要步骤可以分为建立高斯金字塔、生成DOG高斯差分金字塔和DOG局部极值点检测。

　　1.1 尺度空间的构建

　　图像的尺度空间是这幅图像在不同解析度下的表示。一幅图像可以产生几组（octave）图像，一组图像包括几层图像。构造尺度空间传统的方法即构造一个高斯金字塔，原始图像作为最底层，然后对图像进行高斯模糊再降采样（2倍）作为下一层图像（即尺度越大，图像越模糊），循环迭代下去。

对图像进行尺度变换，以满足特征点的尺度不变性，保留图像轮廓和细节。

 

　　1.2 DOG算子

　　DoG（Difference of Gaussian）函数：
　　
　　
　　该函数在计算上只需相邻高斯平滑后图像相减，因此简化了计算。

 

　　1.2.1 DoG高斯差分金字塔

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

 

　　1.3 DOG局部极值检测

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

　　

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

 

2、关键点的精确定位

　　2.1 关键点精确定位

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

 

 　　

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

 　　

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

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

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

 

　　2.2 去除边缘响应

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

 　　

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

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

 　　 

　　该值在两特征值相等时达最小。

 

3、关键点的主方向分配

　　通过尺度不变性求极值点，可以使其具有缩放不变的性质。而利用关键点邻域像素的梯度方向分布特性，可以为每个关键点指定方向参数方向，从而使描述子对图像旋转具有不变性通过求每个极值点的梯度来为极值点赋予方向。
像素点的梯度表示：
　　
梯度幅值：
　　

梯度方向：
　　

 

 

 

 

4、关键点的特征描述

　　下图是一个SIFT描述子事例。其中描述子由2×2×8维向量表征，也即是2×2个8方向的方向直方图组成。左图的种子点由8×8单元组成。每一个小格都代表了特征点邻域所在的尺度空间的一个像素，箭头方向代表了像素梯度方向，箭头长度代表该像素的幅值。然后在4×4的窗口内计算8个方向的梯度方向直方图。绘制每个梯度方向的累加可形成一个种子点，如右图所示：一个特征点由4个种子点的信息所组成。
　　
　　

二、SIFT算法关键点的匹配

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

　　模板图中关键点描述子：

　　实时图中关键点描述子：

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

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

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

　　

 

三、SIFT算法的实现

　　计算图像的SIFT特征，需要用到开源工具包VLFeat。下载地址：www.vlfeat.org。下载后将vl.dll和sift.exe放到项目的目录下。

　　SIFT算法：

from pylab import *
from PIL import Image
import numpy as np
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 = np.loadtxt(filename)
    return f[:, :4], f[:, 4:]  # 特征位置，描述子


def plot_features(im, locs, circle=False):
    """显示带有特征的图像
        输入：im（数组图像），locs（每个特征的行、列、尺度和方向角度）"""

    def draw_circle(c, r):
        t = np.arange(0, 1.01, .01) * 2 * np.pi
        x = r * np.cos(t) + c[0]
        y = r * np.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')
    return


def match(desc1, desc2):
    """对于第一幅图像的每个描述子，选取其在第二幅图像中的匹配
        输入：desc1（第一幅图像中的描述子），desc2（第二幅图像中的描述子）"""

    desc1 = np.array([d / np.linalg.norm(d) for d in desc1])
    desc2 = np.array([d / np.linalg.norm(d) for d in desc2])

    dist_ratio = 0.6
    desc1_size = desc1.shape

    matchscores = np.zeros((desc1_size[0], 1), 'int')
    desc2t = desc2.T    # 预先计算矩阵转置
    for i in range(desc1_size[0]):
        dotprods = np.dot(desc1[i, :], desc2t) # 向量点乘
        dotprods = 0.9999*dotprods
        # 反余弦和反排序，返回第二幅图像中特征的索引
        index = np.argsort(np.arccos(dotprods))

        # 检查最近邻的角度是否小于dist_ratio乘以第二近邻的角度
        if np.arccos(dotprods)[index[0]] < dist_ratio * np.arccos(dotprods)[index[1]]:
            matchscores[i] = int(index[0])

    return matchscores


def match_twosided(desc1, decs2):
    """双向对称版本的match"""

    matches_12 = match(desc1, decs2)
    matches_21 = match(decs2, 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):
    """返回将两幅图像并排拼接成的一幅新图像"""

    # 选取具有最少行数的图像，然后填充足够的空行
    row1 = im1.shape[0]
    row2 = im2.shape[0]

    if row1 < row2:
        im1 = np.concatenate((im1, np.zeros((row2 - row1, im1.shape[1]))), axis=0)
    elif row1 > row2:
        im2 = np.concatenate((im2, np.zeros((row1 - row2, im2.shape[1]))), axis=0)

    # 如果这些情况都没有，那么他们的行数相同，不需要进行填充
    return np.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 = np.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')


if __name__ == '__main__':

    im1f = r'IMG_01.jpg'
    im2f = r'IMG_03.jpg'
    im1 = np.array(Image.open(im1f))
    im2 = np.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()

　　运行结果：

　　

 

 

 

四、SIFT算法的优缺点

4.1 SIFT算法的优点：

　　图像的局部特征，对旋转、尺度缩放、亮度变化保持不变，对视角变化、仿射变换、噪声也保持一定程度的稳定性。

　　独特性好，信息量丰富，适用于海量特征库进行快速、准确的匹配。

　　多量性，即使是很少几个物体也可以产生大量的SIFT特征

　　高速性，经优化的SIFT匹配算法甚至可以达到实时性

　　扩招性，可以很方便的与其他的特征向量进行联合

 

4.2 SIFT算法的缺点：

　　因为是通过对特征点构造128维的向量,然后对向量进行匹配,这样图像就得满足足够多的纹理,否则构造出的128维向量区别性就不是太大,容易造成误匹配,极限情况如指纹图像的匹配,星图识别等这类图像特征点周围根本没有什么纹理,这时SIFT算法就完全失效了。

 

五、地理标记图像匹配

　　由于最后需要对匹配后的图像进行连接可视化，所以需要用到pydot工具包。这个工具包需要先安装配置Graphviz环境才能使用。

　　相关下载地址和安装教程：

　　下载地址：https://graphviz.gitlab.io/_pages/Download/windows/graphviz-2.38.msi

　　教程：https://blog.csdn.net/lanchunhui/article/details/49472949

　　注：最好将图片像素调整到100K左右或者更小，否则可能会出现 xxx.sift not found 或者关于数组维度的报错。像素调小也会让代码运行时间缩短很多。

　　代码实现：

from pylab import *
from PIL import Image
from PCV.localdescriptors import sift
from PCV.tools import imtools
import pydot

# 图像文件夹地址
images_path = "D:/LearnSrc/PythonSrc/Homework/Lecture03/images"
path = "D:/LearnSrc/PythonSrc/Homework/Lecture03/images"

# 获取图像名 形成存放到列表中 计算列表长度
imlist = imtools.get_imlist(images_path)
nbr_images = len(imlist)

# 提取特征值
featlist = [imname[:-3] + 'sift' for imname in imlist]
for i, imname in enumerate(imlist):
    sift.process_image(imname, featlist[i])

matchscores = np.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)

# 复制值
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  # 创建链接所需的最小匹配数

g = pydot.Dot(graph_type='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(pydot.Node(str(i), fontcolor='transparent', shape='rectangle', image=filename))

            # 成对的第二个图像
            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('connect-images.png')

　　运行结果：

 

　　

 

 

内容参考：

https://blog.csdn.net/qq_40369926/article/details/88597406

https://www.pianshen.com/article/1543280501/