SIFT特征提取与检索

1.实验基础

1.1SIFT算法简介

尺度不变特征转换即SIFT (Scale-invariant feature transform)是一种计算机视觉的算法。它用来侦测与描述影像中的局部性特征，它在空间尺度中寻找极值点，并提取出其位置、尺度、旋转不变量，此算法由 David Lowe在1999年所发表，2004年完善总结。
其应用范围包含物体辨识、机器人地图感知与导航、影像缝合、3D模型建立、手势辨识、影像追踪和动作比对。
局部影像特征的描述与侦测可以帮助辨识物体，SIFT特征是基于物体上的一些局部外观的兴趣点而与影像的大小和旋转无关。对于光线、噪声、些微视角改变的容忍度也相当高。基于这些特性，它们是高度显著而且相对容易撷取，在母数庞大的特征数据库中，很容易辨识物体而且鲜有误认。使用 SIFT特征描述对于部分物体遮蔽的侦测率也相当高，甚至只需要3个以上的SIFT物体特征就足以计算出位置与方位。在现今的电脑硬件速度下和小型的特征数据库条件下，辨识速度可接近即时运算。SIFT特征的信息量大，适合在海量数据库中快速准确匹配。
SIFT算法的实质是在不同的尺度空间上查找关键点(特征点)，并计算出关键点的方向。SIFT所查找到的关键点是一些十分突出，不会因光照，仿射变换和噪音等因素而变化的点，如角点、边缘点、暗区的亮点及亮区的暗点等。

1.2算法特点

SIFT算法具有如下特点：

●不变性：SIFT描述图像的局部特征，对图像的旋转和尺度变化具有不变性，对三维视角变化和光照变化具有很强的适应性，在遮挡和场景杂乱时仍保持不变性；

●辨别能力强：特征之间相互区分的能力强，能够在海量的特征中进行快速、准确地匹配；

●特征点数量较多：即使是一幅较小的图片也能够产生足够多的特征点，如一幅500×500的图像能提取出约2000个特征点；

●可扩展性强：可以很方便的与其他形式的特征向量进行联合。

1.3实验步骤

1、尺度空间的极值检测 搜索所有尺度空间上的图像，通过高斯微分函数来识别潜在的对尺度和选择不变的兴趣点。
2、特征点定位 在每个候选的位置上，通过一个拟合精细模型来确定位置尺度，关键点的选取依据他们的稳定程度。
3、特征方向赋值 基于图像局部的梯度方向，分配给每个关键点位置一个或多个方向，后续的所有操作都是对于关键点的方向、尺度和位置进行变换，从而提供这些特征的不变性。
4、特征点描述 在每个特征点周围的邻域内，在选定的尺度上测量图像的局部梯度，这些梯度被变换成一种表示，这种表示允许比较大的局部形状的变形和光照变换。

1.4实验素材

此次实验采用15张图片作为实验素材

2.SIFT特征提取

2.1代码实现

# -*- 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 = '../data/empire.jpg'
im = array(Image.open(imname).convert('L'))
sift.process_image(imname, 'empire.sift')
l1, d1 = sift.read_features_from_file('empire.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()

2.2结果展示

2.3结果分析

1、SIFT特征提取和检测算法提取到的特征点明显多于Harris角点检测，SIFT算法明显优于Harris算法。在一定程度上这体现了SIFT特征提取与检测算法的优越性。
2、从时间效率上来分析，SIFT特征提取和检测消耗了大量的时间，对于一些像素点较多的图像的检测时间可达几小时之久，但其提取特征的准确性不可忽略，大可忽略其时间效率影响。可进行进一步分析，对SIFT算法进行简单优化。

3.SIFT特征匹配

3.1代码实现

 # -*- coding: utf-8 -*-
from pylab import *
from PIL import Image

from PCV.localdescriptors import harris
from PCV.tools.imtools import imresize

"""
This is the Harris point matching example in Figure 2-2.
"""

# Figure 2-2上面的图
#im1 = array(Image.open("../data/crans_1_small.jpg").convert("L"))
#im2= array(Image.open("../data/crans_2_small.jpg").convert("L"))

# Figure 2-2下面的图
im1 = array(Image.open("../data/sf_view1.jpg").convert("L"))
im2 = array(Image.open("../data/sf_view2.jpg").convert("L"))

# resize加快匹配速度
im1 = imresize(im1, (im1.shape[1]/2, im1.shape[0]/2))
im2 = imresize(im2, (im2.shape[1]/2, im2.shape[0]/2))

wid = 5
harrisim = harris.compute_harris_response(im1, 5)
filtered_coords1 = harris.get_harris_points(harrisim, wid+1)
d1 = harris.get_descriptors(im1, filtered_coords1, wid)

harrisim = harris.compute_harris_response(im2, 5)
filtered_coords2 = harris.get_harris_points(harrisim, wid+1)
d2 = harris.get_descriptors(im2, filtered_coords2, wid)

print 'starting matching'
matches = harris.match_twosided(d1, d2)

figure()
gray() 
harris.plot_matches(im1, im2, filtered_coords1, filtered_coords2, matches)
show()

3.2结果展示

3.3结果分析

1、从两租匹配实验可以看出，其图像中的建筑物主体几乎完全匹配，由此可以看出SIFT特征匹配的优越表现。
2、在图片像素不是过高的前提下，SIFT算子匹配的速度很快，甚至可以达到实时的要求。而且可以很方便的与其他形式的特征向量进行联合。

4.特征检索匹配

4.1代码实现

# -*- coding: utf-8 -*-
from PIL import Image
from pylab import *
from numpy import *
import os
from PCV.localdescriptors import sift
from PCV.tools.imtools import get_imlist # 导入原书的PCV模块
import matplotlib.pyplot as plt # plt 用于显示图片
import matplotlib.image as mpimg # mpimg 用于读取图片

# 匹配最多的三张照片

# 获取project2_data文件夹下的图片文件名(包括后缀名)
filelist = get_imlist('project2_data/')

# 输入的图片
im1f = '23.jpg'

im1 = array(Image.open(im1f))
sift.process_image(im1f, 'out_sift_1.txt')
l1, d1 = sift.read_features_from_file('out_sift_1.txt')

i=0
num = [0]*30 # 存放匹配值

for infile in filelist: # 对文件夹下的每张图片进行如下操作

    im2 = array(Image.open(infile))
    sift.process_image(infile, 'out_sift_2.txt')
    l2, d2 = sift.read_features_from_file('out_sift_2.txt')
    matches = sift.match_twosided(d1, d2)
    num[i] = len(matches.nonzero()[0])
    i=i+1
    print '{} matches'.format(num[i-1]) # 输出匹配值
    
i=1
figure()
gray()
while i<4: # 循环三次，输出匹配最多的三张图片

    index=num.index(max(num))
    print index, filelist[index]
    lena = mpimg.imread(filelist[index])  # 读取当前匹配最大值的图片
    # 此时 lena 就已经是一个 np.array 了，可以对它进行任意处理
    subplot(1,3,i)
    plt.imshow(lena)  # 显示图片
    plt.axis('off')  # 不显示坐标轴
    num[index] = 0  #将当前最大值清零
    i=i+1
    
show()

4.2结果展示

输入

输出

4.3结果分析

1、输入图片与数据集中图片两两进行特征匹配后得到的匹配值保存在数组num中，再利用max()函数以及index()函数分别得到匹配最大值和匹配最大值的下标，之后利用下标，通过matplotlib库即可显示照片，要找到与输入图片匹配特征点最多的图，最后将矩阵中数值最大的三个位置输出，这里用了一个while循环，在每次循环结束时将当前循环所得到的匹配最大值清零，这样才不会得到重复的答案。SIFT在大量数据库中速度也非常快。
2、有实验结果可以看出，输出结果基本上按照相似度由高到低排序，数据匹配结果精确，说明sift算法在特征明显的复杂场景下运用效果更好。

5.匹配地理标记

5.1代码实现

# -*- 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:\Pictures\one"  # set this to the path where you downloaded the panoramio images
path = "F:\Pictures\one"  # 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('11.jpg')

5.2结果展示

5.3结果分析

1、在涵盖较多复杂场景的15张图片集，通过算法匹配，进而排序，可得观察到其得出的实验结果与我们目视判读的结果一致，一次可以得出SIFT特征匹配的准确率是相当高的。
2、该实验耗时较长，大一些的图片或者是复杂一点的图片运行时间越长，建议缩小图片像素大小以增强实时性。
3、实验效果好与高效提取关键点、准确地对关键点附加详细的信息(局部特征)、通过特征点(附带上特征向量的关 键点)的两两比较找出相互匹配的若干对特征点,建立景物间的对应关系有密不可分地关系。

6.RANSAC算法

RANSAC基本原理

RANSAC为Random Sample Consensus的缩写，它是根据一组包含异常数据的样本数据集，计算出数据的数学模型参数，得到有效样本数据的算法。它于1981年由Fischler和Bolles最先提出 。该算法在计算机视觉中广泛应用。

RANSAC算法的基本假设是样本中包含正确数据(inliers，可以被模型描述的数据)，也包含异常数据(outliers，偏离正常范围很远、无法适应数学模型的数据)，即数据集中含有噪声。这些异常数据可能是由于错误的测量、错误的假设、错误的计算等产生的。同时RANSAC也假设，给定一组正确的数据，存在可以计算出符合这些数据的模型参数的方法。 因此，RANSAC通过反复选择数据中的一组随机子集来达成目标。被选取的子集被假设为局内点，并用下述方法进行验证：

基本RANSAC算法流程

1、随机抽取部分数据认为选中的数据是inliers, 并构建数学模型。
2、用1中得到的模型去测试所有的其它数据，如果某个点适用于估计的模型，认为它也是局内点，否是就是局外点。
3、统计1中模型的局内点数目。
4、重复1~3给定的次数，维护全局最值，认为局内点最多的模型为我们最终需要的模型。

RANSAC与最小二乘法的区别

最小二乘法尽量去适应包括局外点在内的所有点。相反，RANSAC能得出一个仅仅用局内点计算出模型，并且概率还足够高。但是，RANSAC并不能保证结果一定正确，为了保证算法有足够高的合理概率，必须小心的选择算法的参数（参数配置）。经实验验证，对于包含80%误差的数据集，RANSAC的效果远优于直接的最小二乘法。

实现代码

# -*- coding: utf-8 -*-
import cv2
import numpy as np
import random

def compute_fundamental(x1, x2):
    n = x1.shape[1]
    if x2.shape[1] != n:
        raise ValueError("Number of points don't match.")

    # build matrix for equations
    A = np.zeros((n, 9))
    for i in range(n):
        A[i] = [x1[0, i] * x2[0, i], x1[0, i] * x2[1, i], x1[0, i] * x2[2, i],
                x1[1, i] * x2[0, i], x1[1, i] * x2[1, i], x1[1, i] * x2[2, i],
                x1[2, i] * x2[0, i], x1[2, i] * x2[1, i], x1[2, i] * x2[2, i]]

    # compute linear least square solution
    U, S, V = np.linalg.svd(A)
    F = V[-1].reshape(3, 3)

    # constrain F
    # make rank 2 by zeroing out last singular value
    U, S, V = np.linalg.svd(F)
    S[2] = 0
    F = np.dot(U, np.dot(np.diag(S), V))

    return F / F[2, 2]


def compute_fundamental_normalized(x1, x2):
    """    Computes the fundamental matrix from corresponding points
        (x1,x2 3*n arrays) using the normalized 8 point algorithm. """

    n = x1.shape[1]
    if x2.shape[1] != n:
        raise ValueError("Number of points don't match.")

    # normalize image coordinates
    x1 = x1 / x1[2]
    mean_1 = np.mean(x1[:2], axis=1)
    S1 = np.sqrt(2) / np.std(x1[:2])
    T1 = np.array([[S1, 0, -S1 * mean_1[0]], [0, S1, -S1 * mean_1[1]], [0, 0, 1]])
    x1 = np.dot(T1, x1)

    x2 = x2 / x2[2]
    mean_2 = np.mean(x2[:2], axis=1)
    S2 = np.sqrt(2) / np.std(x2[:2])
    T2 = np.array([[S2, 0, -S2 * mean_2[0]], [0, S2, -S2 * mean_2[1]], [0, 0, 1]])
    x2 = np.dot(T2, x2)

    # compute F with the normalized coordinates
    F = compute_fundamental(x1, x2)
    # print (F)
    # reverse normalization
    F = np.dot(T1.T, np.dot(F, T2))

    return F / F[2, 2]

def randSeed(good, num = 8):
    '''
    :param good: 初始的匹配点对
    :param num: 选择随机选取的点对数量
    :return: 8个点对list
    '''
    eight_point = random.sample(good, num)
    return eight_point

def PointCoordinates(eight_points, keypoints1, keypoints2):
    '''
    :param eight_points: 随机八点
    :param keypoints1: 点坐标
    :param keypoints2: 点坐标
    :return:8个点
    '''
    x1 = []
    x2 = []
    tuple_dim = (1.,)
    for i in eight_points:
        tuple_x1 = keypoints1[i[0].queryIdx].pt + tuple_dim
        tuple_x2 = keypoints2[i[0].trainIdx].pt + tuple_dim
        x1.append(tuple_x1)
        x2.append(tuple_x2)
    return np.array(x1, dtype=float), np.array(x2, dtype=float)


def ransac(good, keypoints1, keypoints2, confidence,iter_num):
    Max_num = 0
    good_F = np.zeros([3,3])
    inlier_points = []
    for i in range(iter_num):
        eight_points = randSeed(good)
        x1,x2 = PointCoordinates(eight_points, keypoints1, keypoints2)
        F = compute_fundamental_normalized(x1.T, x2.T)
        num, ransac_good = inlier(F, good, keypoints1, keypoints2, confidence)
        if num > Max_num:
            Max_num = num
            good_F = F
            inlier_points = ransac_good
    print(Max_num, good_F)
    return Max_num, good_F, inlier_points


def computeReprojError(x1, x2, F):
    """
    计算投影误差
    """
    ww = 1.0/(F[2,0]*x1[0]+F[2,1]*x1[1]+F[2,2])
    dx = (F[0,0]*x1[0]+F[0,1]*x1[1]+F[0,2])*ww - x2[0]
    dy = (F[1,0]*x1[0]+F[1,1]*x1[1]+F[1,2])*ww - x2[1]
    return dx*dx + dy*dy

def inlier(F,good, keypoints1,keypoints2,confidence):
    num = 0
    ransac_good = []
    x1, x2 = PointCoordinates(good, keypoints1, keypoints2)
    for i in range(len(x2)):
        line = F.dot(x1[i].T)
        #在对极几何中极线表达式为[A B C],Ax+By+C=0,  方向向量可以表示为[-B,A]
        line_v = np.array([-line[1], line[0]])
        err = h = np.linalg.norm(np.cross(x2[i,:2], line_v)/np.linalg.norm(line_v))
        # err = computeReprojError(x1[i], x2[i], F)
        if abs(err) < confidence:
            ransac_good.append(good[i])
            num += 1
    return num, ransac_good


if __name__ =='__main__':
    im1 = r'F:\Pictures\one\3.jpg'
    im2 = r'F:\Pictures\one\03.jpg'

    print(cv2.__version__)
    psd_img_1 = cv2.imread(im1, cv2.IMREAD_COLOR)
    psd_img_2 = cv2.imread(im2, cv2.IMREAD_COLOR)
    # 3) SIFT特征计算
    sift = cv2.xfeatures2d.SIFT_create()
    # find the keypoints and descriptors with SIFT
    kp1, des1 = sift.detectAndCompute(psd_img_1, None)
    kp2, des2 = sift.detectAndCompute(psd_img_2, None)

    # FLANN 参数设计
    match = cv2.BFMatcher()
    matches = match.knnMatch(des1, des2, k=2)

    # Apply ratio test
    # 比值测试，首先获取与 A距离最近的点 B （最近）和 C （次近），
    # 只有当 B/C 小于阀值时（0.75）才被认为是匹配，
    # 因为假设匹配是一一对应的，真正的匹配的理想距离为0
    good = []
    for m, n in matches:
        if m.distance < 0.75 * n.distance:
            good.append([m])
    print(good[0][0])

    print("number of feature points:",len(kp1), len(kp2))
    print(type(kp1[good[0][0].queryIdx].pt))
    print("good match num:{} good match points:".format(len(good)))
    for i in good:
        print(i[0].queryIdx, i[0].trainIdx)


    Max_num, good_F, inlier_points = ransac(good, kp1, kp2, confidence=30, iter_num=500)
    # cv2.drawMatchesKnn expects list of lists as matches.
    # img3 = np.ndarray([2, 2])
    # img3 = cv2.drawMatchesKnn(img1, kp1, img2, kp2, good[:10], img3, flags=2)

    # cv2.drawMatchesKnn expects list of lists as matches.

    img3 = cv2.drawMatchesKnn(psd_img_1,kp1,psd_img_2,kp2,good,None,flags=2)
    img4 = cv2.drawMatchesKnn(psd_img_1,kp1,psd_img_2,kp2,inlier_points,None,flags=2)
    cv2.namedWindow('image1', cv2.WINDOW_NORMAL)
    cv2.namedWindow('image2', cv2.WINDOW_NORMAL)
    cv2.imshow("image1",img3)
    cv2.imshow("image2",img4)
    cv2.waitKey(0)#等待按键按下
    cv2.destroyAllWindows()#清除所有窗口

对比试验

（1）实验一：景深单一

sift特征提取匹配点：

利用RANSAC算法剔除错误匹配点：

小结：通过对比实验可以发现，SIFT特征匹配出现了很多错误匹配点，而RAMSAC算法去除了一部分干扰项，剔除错误匹配点的效果很好。在景深单一，特征点比较明显的情况下，RANSAC算法剔除错误匹配点的效果很好。因为图像没有旋转和太大的角度变换， SIFT特征匹配能够匹配出较多的匹配点，但是利用RANSAC算法之后，可以明显看到有的匹配正确的点也会被删除掉，留下的少部分的匹配点完全正确。 所以经过筛选后的匹配点虽然位置很精准，但是同时也减弱了SIFT特征匹配的效果。

（2）实验二：景深丰富

sift特征提取匹配点：

利用RANSAC算法剔除错误匹配点：

小结：由该组对比试验可以看出，在景深丰富的情况下，SIFT特征匹配在近景产生的匹配点较多，且有一些错误匹配点，相比较之下远景产生的匹配点较少。进行RANSAC算法剔除之后，留下了的匹配点比较精准。和实验一一样，RANSAC算法删除了一些正确的匹配点，剩下的匹配点数量较少

7.实验总结

1、与之前所学的Harris算法相比，SIFT特征不只具有尺度不变性，即使改变旋转角度，图像亮度或拍摄视角，仍然能够得到好的检测效果，这是相比Harris算法更优的地方。
2、sift算法较为复杂，使用像素点过高的图片会导致运行时间过长，因此建议使用像素点少的图片以达到更快的实验结果。
3、SIFT算法独特性好，信息量丰富，适用于在海量特征数据库中进行快速、准确的匹配。在实验过程中有明确的体现，尤其在SIFT算法与Harris算法的对比中可以看出SIFT算法检测到的信息量丰富。
4、实验过程中遇到的问题：
（1）运行速度过长：由于图像像素过大导致运行迟迟未结束，经过对图像进行压缩，压缩后的图片再运行可以显著提高速度。
（2）无法进行特征匹配：由于图片经过了压缩裁剪，所以尺寸有些许不同，可是该算法只能对尺寸相同图片进行匹配，因此需要使用PS将数据集中图片改成相同大小尺寸。