OpenCV+Python实现SIFT匹配地理标记图像

本文将持续更新 21.3.29

1. 什么是SIFT

SIFT的全称是Scale Invariant Feature Transform，尺度不变特征变换，2004年，不列颠哥伦比亚大学的D.Lowe在他的论文“Distinctive Image Features from Scale-Invariant Keypoints”（尺度不变关键点中的独特图像特征）中提出了一种新的尺度不变特征变换（SIFT）算法，该算法提取关键点并计算其描述符。
SIFT特征对旋转、尺度缩放、亮度变化等保持不变性，是一种非常稳定的局部特征。

1.1 SIFT算法解决的问题

以往的角点检测算法，例如Harris角点检测器可以做到旋转不变（rotation-invariant），即无论图像发生旋转，算法依然可以找到相同的角。但是当图像发生缩放之后，原先的角可能就发生了变化，算法就无法检测出角点。例如下图，在同一检测窗口中放大小窗口中小图像中的角时，该角是平坦的。为了解决这个问题，我们需要一种尺度不变（或者叫缩放不变）算法。

1.2 SIFT算法的步骤

1.2.1 尺度空间极值检测
如上图，对于一个图像我们不能使用同样大小的检测窗口来检测不同放缩比例的关键点。对一个“小”角点来说是没问题啦，但要检测更“大”的角点，我们需要更大的窗口。为此，使用尺度空间过滤。

在其中，不同的σ值的图像对应的高斯-拉普拉斯算子被算出来。LoG（Laplacian of Gaussian，详见高斯拉普拉斯算子-知乎就是前文的高斯-拉普拉斯算子）作为一个斑点检测器，能检测出由于σ变化引起的不同大小下的斑点。

简单的说，σ 就是一个放缩比例的参数。例如，在上图中，有较小 σ 值的高斯内核对一个“小”角点返回一个较大的值，然而较大的 σ值的高斯内核适合较“大”的角点。所以，我们可以在放缩空间中求出局部最大值，它给我们返回一个(x,y,σ)值的列表，它告诉我们，在(x,y)点 σ 放缩比的情况下，有一个潜在的关键点。
原文：https://blog.csdn.net/ssybc/article/details/86628410

但是LoG比较消耗性能，因此 SIFT 算法使用了高斯差分算子(Difference of Gaussian，高斯差分算子)，它是LoG的近似值。

高斯差分算子可以通过具有两个不同的图像的高斯模糊的差来获得σ。令他们为 σ 和 kσ，这个过程为不同图像的高斯金字塔的阶进行了处理。
如下图所示：

1.2.2 高斯尺度空间

我们还可以通过图像的模糊程度来模拟人在距离物体由远到近时物体在视网膜上成像过程，距离物体越近其尺寸越大图像也越模糊，这就是高斯尺度空间，使用不同的参数模糊图像（分辨率不变），是尺度空间的另一种表现形式。
我们知道图像和高斯函数进行卷积运算能够对图像进行模糊，使用不同的“高斯核”可得到不同模糊程度的图像。一副图像其高斯尺度空间可由其和不同的高斯卷积得到：
 
L(x,y,σ)=G(x,y,σ)∗I(x,y)L(x,y,σ)=G(x,y,σ)∗I(x,y)
 
其中，G(x,y,σ)是高斯核函数。G(x,y,σ)是高斯核函数。

 
G(x,y,σ)=12πσ2ex2+y22σ2G(x,y,σ)=12πσ2ex2+y22σ2
 

σσ称为尺度空间因子，它是高斯正态分布的标准差，反映了图像被模糊的程度，其值越大图像越模糊，对应的尺度也就越大。L(x,y,σ)L(x,y,σ)代表着图像的高斯尺度空间。
构建尺度空间的目的是为了检测出在不同的尺度下都存在的特征点，而检测特征点较好的算子是Δ2GΔ2G(高斯拉普拉斯,LoG）

 
Δ2=∂2∂x2+∂2∂y2Δ2=∂2∂x2+∂2∂y2
 

1.2.3 高斯金字塔构建示例
以一个512×512512×512的图像I为例，构建高斯金字塔步骤：(从0开始计数，倒立的金字塔）

金字塔的组数，log2512=9log2⁡512=9，减去因子3，构建的金字塔的组数为6。取每组的层数为3。
构建第0组，将图像的宽和高都增加一倍，变成1024×10241024×1024（I0I0）。第0层I0∗G(x,y,σ0)I0∗G(x,y,σ0)，第1层I0∗G(x,y,kσ0)I0∗G(x,y,kσ0)，第2层I0∗G(x,y,k2σ0)I0∗G(x,y,k2σ0)
构建第1组，对I0I0降采样变成512×512512×512（I1I1）。第0层I1∗G(x,y,2σ0)I1∗G(x,y,2σ0)，第1层I1∗G(x,y,2kσ0)I1∗G(x,y,2kσ0)I1∗G(x,y,2k2σ0)I1∗G(x,y,2k2σ0)
⋮⋮
构建第o组，第s层 Io∗G(x,y,2oksσ0)Io∗G(x,y,2oksσ0)
高斯金字塔构建成功后，将每一组相邻的两层相减就可以得到DoG金字塔.

算出了 DoG 之后，图像就开始在放缩空间中搜索局部极值。比方说，图像中的某像素点用来和它周围的8个邻居比较。
例如，将一幅图像中的一个像素与它的8个相邻像素以及下一个比例中的9个像素和前一个比例中的9个像素进行比较，如果它是局部的极大值，它就是一个潜在的关键点。基本上意思就是说那个关键点，在此比例下是最有代表性的。
如下图所示：

对于不同的参数，论文给出了一些经验数据，可以概括为，阶数=4，放缩等级数=5，初始σ=1.6，k=\sqrt{2}等作为最优值。
原文：https://blog.csdn.net/ssybc/article/details/86628410

2.关键点定位
一旦潜在关键点的位置被算出来了，为了得到更精确的结果，它们必须被再精选一次。他们使用放缩空间中的泰勒级数展开式，得到了更精确的极值点位置。并且如果该极值处的强度小于阈值(该论文为0.03)，则不判定为关键点。

这个阈值在OpenCV里叫做contrastThreshold。

DoG 算法对于边缘也有很高的响应，所以我们得把边缘除去。对此，使用了一个与哈里斯角点检测相似的概念。他们使用一个 2x2 的黑塞矩阵(H)(译者附：黑塞矩阵)来计算主曲率。由Harris角检测器可知，对于边，一个特征值远大于另一个特征值。

这里用了一个简单的函数，如果该比例大于在OpenCV里被称为edgeThreshold的阈值，该关键点会被抛弃。在论文中该阈值被设为10。

因此，它消除了任何低对比度的关键点和边缘关键点，剩下的是强烈的兴趣点。

3. 方向指定
现在为每个关键点指定一个方向，以实现图像旋转的不变性。邻域取围绕关键点位置的区域，大小取决于放缩比例。在那个区域计算出梯度大小和梯度方向。
一个有着36个抽屉、覆盖360度的方向直方图被创建了出来(它通过梯度幅值和高斯加权圆形窗口来加上了权重，其中σ取关键点的1.5倍放缩比)。取该直方图中最高的峰值，以及任何在最高峰值80%之上的峰值点，都用于计算方向。它创建的关键点具有相同的位置和放缩比，但方向不同。有利于匹配的稳定性。

4. 关键点描述
现在关键点的描述符已经创建好了。取围绕关键点的一个 16x16 邻域，它被分成16个4x4大小的子块。对于每个子块，创建8抽屉的方向直方图，这样总共就有了128个可以用的有值得抽屉。它用一个向量来表示，形成了关键点描述符。除此之外，还采取了一些措施来提升对抗光照变化、旋转等的健壮性。

5. 关键点匹配
两张图像之间的关键点，通过识别它们之间的最近的邻居来匹配。但在某些情况之下，第二相近的匹配度有可能和最相近的匹配度非常接近。这有可能是因为噪音或者其他某些原因引起的。如果是那样的话，取次接近的和最接近的比值。如果比值大于0.8，它们就被丢弃掉。在该论文中使用的这种方案，丢弃掉了大约 90% 的错误匹配，但只丢掉了5%的正确匹配。

这就是SIFT算法的一个简单的摘要。想要理解更多的细节，强烈推荐你阅读原论文。

2. SIFT实际应用

2.1 OpenCV中的SIFT

这是根据Opencv文档中的示例写的一个演示代码，和原文不同的地方就是SIFT调用方式：“sift = cv.SIFT_create()” 需要改为 “sift = cv.xfeatures2d.SIFT_create()”

import cv2 as cv

img = cv.imread(r'test_pic\soc.jpg',1)
sift = cv.xfeatures2d.SIFT_create()
kp = sift.detect(img,None)
img=cv.drawKeypoints(img,kp,img)
cv.imwrite(r'test_pic\sift_keypoints.jpg',img)

结果如下：

2.SIFT特征匹配

好了，你现在已经知道了SIFT的基本原理，开始实现SIFT特征匹配吧！（笑）

实现特征匹配我们需要通过SIFT查找得到关键点和特征向量，

import cv2 as cv
import matplotlib.pyplot as plt

img1 = cv.imread(r'test_pic\cool1.jpg')
img2 = cv.imread(r'test_pic\cool2.jpg')
# 创建SIFT特征点检测
sift = cv.xfeatures2d.SIFT_create()
# 计算关键点和特征符
kp1, des1 = sift.detectAndCompute(img1, None)
kp2, des2 = sift.detectAndCompute(img2, None)

# BFmatcher with default parms
bf = cv.BFMatcher(cv.NORM_L2)
matches = bf.knnMatch(des1, des2, k=2) 
# 储存特征向量匹配最好的优质匹配点
goodMatchs = []
for m, n in matches:
    if m.distance < 0.7 * n.distance:
        goodMatchs.append(m)
        
pic3 = cv.drawMatches(img1=img1, keypoints1=kp1, img2=img2, keypoints2=kp2, matches1to2=goodMatchs, outImg=None)
cv.imwrite(r'test_pic\compare.jpg',pic3)

结果如下：

3.SIFT算法匹配地理标记图像

1.匹配地理标记图像

2.代码

1.前期准备一下缩略图
将缩略图的保存和地理标记图像匹配的主体功能分开的原因是：OpenCV的imwrite()函数会在py文件运行结束后，才将图片保存到指定路径，这将导致pydotplus和graphviz画出的无向图找不到缩略图而报错:

save_pic.py：

import os
import cv2 as cv
from pylab import *

maxsize = (100, 100)  # 定义缩略图的大小
path = r'E:\PycharmProjects\SIFT\pic\result\pic'

# 读取整个文件夹的图片
def read_path(pathname):
    imgname_list = os.listdir(pathname)
    img_list = []
    i = 0
    # 图片列表
    for imgname in imgname_list:
        if imgname.endswith('.jpg'):
            img = cv.imread(pathname + '/' + imgname)
            img_n = cv.resize(img, maxsize, cv.INTER_AREA)
            filename = path + str(i) + '.png'
            cv.imwrite(filename, img_n)  # need temporary files of the right size
            i = i + 1
            print(i)
    return img_list

list = read_path(r'E:\PycharmProjects\SIFT\pic')
print(list)

2.主要的代码如下
SIFT.py：

import os
import cv2 as cv
from pylab import *
import pydotplus as pydot

# 读取整个文件夹的图片
def read_path(pathname):
    imgname_list = os.listdir(pathname)
    img_list = []
    # 图片列表
    for imgname in imgname_list:
        if imgname.endswith('.jpg'):
            img = cv.imread(pathname + '/' + imgname)
            img_list.append(img)
    return img_list


img_list = read_path(r'E:\PycharmProjects\SIFT\pic')
nbr_images = len(img_list)
match_scores = zeros((nbr_images, nbr_images))

for i in range(nbr_images):
    for j in range(i, nbr_images):  # only compute upper triangle
        print('comparing ', i, j)
        sift = cv.xfeatures2d.SIFT_create()
        kp1, des1 = sift.detectAndCompute(img_list[i], None)
        kp2, des2 = sift.detectAndCompute(img_list[j], None)
        # BFMatch匹配
        bf = cv.BFMatcher(cv.NORM_L2)
        matches = bf.knnMatch(des1, des2, k=2)
        # 储存差距小的优秀匹配点
        goodMatches = []
        for m, n in matches:
            if m.distance < 0.5 * n.distance:
                goodMatches.append(m)
        # 计算优秀匹配点的和
        nbr_matches = len(goodMatches)
        # 向match_scores赋值
        print('number of matches = ', nbr_matches)
        match_scores[i, j] = nbr_matches

# 复制
for i in range(nbr_images):
    for j in range(i + 1, nbr_images):  # no need to copy diagonal 不用复制自我匹配的对角线
        match_scores[j, i] = match_scores[i, j]
# 可视化
threshold = 2  #  至少2个以上匹配点就可以算是联系
g = pydot.Dot(graph_type='graph')  # 不需要有向图
maxsize = (100, 100)  # 定义缩略图的大小
path = r'E:\PycharmProjects\SIFT\pic\result\pic'
#两两配对
for i in range(nbr_images):
    for j in range(i + 1, nbr_images):
        if match_scores[i, j] > threshold:
            filename = path + str(i) + '.png'
            g.add_node(pydot.Node(str(i), fontcolor='transparent', shape='rectangle', image=filename))
            filename = path + str(j) + '.png'
            g.add_node(pydot.Node(str(j), fontcolor='transparent', shape='rectangle', image=filename))
            g.add_edge(pydot.Edge(str(i), str(j)))
#绘制S地理标记SIFT匹配图
g.write_jpg('sift.jpg')

3.分步解释
SIFT中每个特征点的描述子都有几个属性：坐标、尺度和梯度方向角度和128维的邻域子向量。

为什么m.distance设定在小于等于0.5？
m.distance是sift特征点描述中128维邻域子向量的归一化欧几里得距离（欧氏距离），根据SIFT的创造者Lowe在论文中提出来的稳健匹配准则，我们应当使用两个特征距离和两个最匹配特征距离的比率。

当我们得到这个欧式距离的比率之后，应该将它们按升序排列，这意味着图像1中的关键点与图像2中的关键点之间的距离是最近的。
之后我们就应该设置一个阈值T，如果m.distance和n.distance的比率小于这个阈值T，那么它就是一个好的匹配，我们将它记录到good_matches里面。

所以这个阈值T，应该被设置的小一点，一般取0.3~0.7。