计算机视觉二（python）——SIFT地理特征匹配

SIFT地理特征匹配

一 SIFT算法：
SIFT(Scale-Invariant Feature Transform,尺度不变特征变换)，由David G.Lowe提出，是一种计算机视觉的算法。它用来侦测与描述影像中的局部性特征，它在空间尺度中寻找极值点，并提取出其位置、尺度、旋转不变量。SIFT算法可解决的问题有目标的旋转、缩放、平移（RST）、图像仿射/投影变换（视点viewpoint）、弱光照影响（illumination）、部分目标遮挡（occlusion）、杂物场景（clutter）、噪声等。

1.1 SIFT特征检测的步骤

（1）尺度空间的极值检测：搜索所有尺度空间上的图像，通过高斯微分函数来识别潜在的对尺度和选择不变的兴趣点。
（2）特征点定位 在每个候选的位置上，通过一个拟合精细模型来确定位置尺度，关键点的选取依据他们的稳定程度。
（3）特征方向赋值 基于图像局部的梯度方向，分配给每个关键点位置一个或多个方向，后续的所有操作都是对于关键点的方向、尺度和位置进行变换，从而提供这些特征的不变性。
（4）特征点描述 在每个特征点周围的邻域内，在选定的尺度上测量图像的局部梯度，这些梯度被变换成一种表示，这种表示允许比较大的局部形状的变形和光照变换。
1.2 SIFT算法原理
sift算法的实质是在不同的尺度空间上查找关键点(特征点)，并计算出关键点的方向。这些是一些十分突出的点，不会因光照、尺度、旋转等因素的改变而消失，比如角点、边缘点、暗区域的亮点以及亮区域的暗点。

1.3 高斯模糊
SIFT算法是在不同的尺度空间上查找关键点，而尺度空间的获取需要使用高斯模糊来实现。高斯核是唯一可以产生多尺度空间的核，一个图像的尺度空间，L(x, y, σ) ,定义为原始图像 I(x, y)与一个可变尺度的二维高斯函数G(x, y, σ) 卷积运算。

在二维空间中，这个公式生成的曲面的等高线是从中心开始呈正态分布的同心圆，如下图所示。分布不为零的像素组成的卷积矩阵与原始图像做变换。每个像素的值都是周围相邻像素值的加权平均。原始像素的值有最大的高斯分布值，所以有最大的权重，相邻像素随着距离原始像素越来越远，其权重也越来越小。这样进行模糊处理比其它的均衡模糊滤波器更高地保留了边缘效果。

1.4 尺度空间
尺度空间主要思想是通过对原始图像进行尺度变换，获得图像多尺度下的空间表示。从而实现边缘、角点检测和不同分辨率上的特征提取，以满足特征点的尺度不变性。
1.5 构建高斯金字塔
图像的金字塔模型是指，将原始图像不断降阶采样，得到一系列大小不一的图像，由大到小，从下到上构成的塔状模型。原图像为金子塔的第一层，每次降采样所得到的新图像为金字塔的一层(每层一张图像)，每个金字塔共n层。金字塔的层数根据图像的原始大小和塔顶图像的大小共同决定。

二、实现地理图像匹配

1.Harris角点检测
Harris角点检测基本思想是从图像局部的小窗口观察图像特征（角点：窗口向任意方向的移动都导致图像灰度的明显变化）。
将图像转换为灰度图像，计算响应函数，基于响应值选取角点，最后在原始图像中覆盖绘制检测出的角点。
Harris角点检测代码如下：

在这from matplotlib import pyplot as plt
import cv2
import numpy as np


def pltshow(imgs, titles):
    for i in range(len(imgs)):
        plt.subplot(2, len(imgs) / 2 + 1, i + 1), plt.imshow(imgs[i], 'gray')
        plt.title(titles[i])
        plt.xticks([]), plt.yticks([])
    plt.show()


def cvshow(name, img):
    cv2.imshow(name, img)
    cv2.waitKey(0)
    cv2.destroyAllWindows()


def getimg(path, flag):
    img = cv2.imread(path)
    b, g, r = cv2.split(img)
    img_c = cv2.merge([r, g, b])
    img_gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
    if flag == 0:
        return img
    # 若要用Matplot显示就先把bgr转换成rgb,img_c显示即可，opencv imshow时仍用img
    elif flag == 1:
        return img_c
    elif flag == 2:
        return img_g


# 图像特征harris角点检测
img = getimg('D:/Software/Pictures/ludaa.jpg', 0)
print('img.shape', img.shape)
gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
dst = cv2.cornerHarris(gray, 2, 3, 0.04)
print('dst.shape', dst.shape)

# 做一个非极大值抑制，
# 对比几种情况，当dst越大的时候，角点更容易出现，否则边界会出现很多
orig1 = img.copy()
orig2 = img.copy()
orig3 = img.copy()
orig1[dst > 0.02 * dst.max()] = [0, 0, 255]
orig2[dst > 0.1 * dst.max()] = [0, 0, 255]
orig3[dst > 0.2 * dst.max()] = [0, 0, 255]

res = np.hstack((orig1, orig2, orig3))
cvshow('dst', res)里插入代码片

运行结果如下：

结果分析：
通过运行结果截图可以看出来，角点的周围还会存在其他的角点，因为角点周围会有响应值比较大的点所以角点的分布会比较密集。

SIFT特征匹配

代码如下：

import numpy as np
import cv2
from matplotlib import pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

img1 = cv2.imread(
    "D:/Software/Pictures/ludae.jpg", 0)  # 导入灰度图像
img2 = cv2.imread(
    "D:/Software/Pictures/ludab.jpg", 0)

detector = cv2.ORB_create()  # Oriented FAST and Rotated BRIEF

# detector = cv2.AKAZE_create()

kp1 = detector.detect(img1, None)
kp2 = detector.detect(img2, None)
kp1, des1 = detector.compute(img1, kp1) # keypoint 是关键点的列表,desc 检测到的特征的局部图的列表
kp2, des2 = detector.compute(img2, kp2)
# 获得一个暴力匹配器的对象
bf = cv2.BFMatcher(cv2.NORM_HAMMING, crossCheck=True)


# 利用匹配器 匹配两个描述符的相近程度
"""knn法"""
matches = bf.knnMatch(des1, des2, k=1)  # knn 匹配可以返回k个最佳的匹配项、bf返回所有的匹配项
img32 = cv2.drawMatchesKnn(img1, kp1, img2, kp2, matches[: 50], img2, flags=2)

plt.figure(figsize=(10,10))
plt.title('knn法',fontsize=12,color='r')
plt.imshow(img32)
plt.axis('off')
plt.show()

cv2.waitKey(0)
cv2.destroyAllWindows()

结果如下：

匹配地理标记图像

代码如下：

```python
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 = 'C:/Users/Mark/project1/task/test'  # set this to the path where you downloaded the panoramio images
path = 'C:/Users/Mark/project1/task/test/'  # 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 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 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('jmu.png')

结果如下：