SIFT特征提取与检索

目录

1.前言

2.兴趣点

3.描述子

4.关键点匹配

5.本次实验图片集如下

6.sift特征点

7.SIFT特征匹配

8.Harris特征匹配

9.多张图片索引出匹配度最高的

10.实验分析

11.遇到的问题与解决方法

---

1.前言

所谓局部特征描述算子（SIFT）就是用来刻画图像中的这些局部共性的，而我们也可以将一幅图像映射（变换）为一个局部特征的集合。理想的局部特征应具有平移、缩放、旋转不变性，同时对光照变化、仿射及投影影响也应有很好的鲁棒性。传统的局部特征往往是直接提取角点或边缘，对环境的适应能力较差。

SIFT针对局部特征进行特征提取，在尺度空间寻找极值点，提取位置，尺度，旋转不变量，生成特征描述子。SIFT算法的实质可以归为在不同尺度空间上查找关键点（特征点）的问题。所谓关键点，就是一些十分突出的点，这些点不会因光照条件的改变而消失，比如角点、边缘点、暗区域的亮点以及亮区域的暗点，既然两幅图像中有相同的景物， 那么使用某种方法分别提取各自的稳定点，这些点之间就会有相互对应的匹配点。而在SIFT中，关键点是在不同尺度空间的图像下检测出的具有方向信息的局部极值点。

 

2.兴趣点

SIFT 特征使用高斯差分函数来定位兴趣点：
                                

其中，  是上一章中介绍的二维高斯核， 是使用  模糊的灰度图像， κ 是决定相差尺度的常数。兴趣点是在图像位置和尺度变化下 D(x,σ) 的最大值和最小值点。这些候选位置点通过滤波去除不稳定点。基于一些准则，比如认为低对比度和位于边上的
点不是兴趣点，我们可以去除一些候选兴趣点。

3.描述子

为了对图像亮度具有稳健性，SIFT 描述子使用图像梯度。 SIFT 描述子在每个像素点附近选取子区域网格，在每个子区域内计算图像梯度方向直方图。每个子区域的直方图拼接起来组成描述子向量。 SIFT 描述子的标准设置使用 4× 4 的子区域，每个子区域使用 8 个小区间的方向直方图，会产生共128 个小区间的直方图（4× 4× 8=128）。如下图所示

4.关键点匹配

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

模板图中关键点描述子：

实时图中关键点描述子：

任何两描述子相似性度量：
要得到配对的关键点描述子，需要满足：

 

关键点的匹配一般都采用kd树（Kd树是一个平衡二叉树）的数据结构来完成搜索。搜索的内容是以目标图像的关键点为基准，搜索与目标图像的特征点最邻近的原图像特征点和次邻近的原图像特征点。

5.本次实验图片集如下

 

6.sift特征点

（1）此处将使用开源工具包VLFeat提供的二进制文件来计算图像的SIFT特征。VLFeat工具包可以从http://www.vlfeat.org/下载，二进制文件可以在所有主要的平台上运行。VLFeat库是用C语言来写的，但是我们可以使用该库提供的命令行接口。以在Windows 10 64bit平台上为例，下载的文件为vlfeat-0.9.20-bin.tar.gz，解压缩后，将vlfeat-0.9.20/bin/win64文件夹下的vl.dll和sift.exe和vl.lib复制粘贴到项目目录下，这样基本配置完成。

（2）sift.py代码如下：

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

def write_featrues_to_file(filename, locs, desc):
    """将特征位置和描述子保存到文件中"""
    savetxt(filename, hstack((locs,desc)))

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

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

imname = 'D:\JMU\computer_vision\SIFTcode\SIFT\scene1.jpg'
im1 = array(Image.open(imname).convert('L'))
process_image(imname, '1.sift')
l1,d1 = read_features_from_file('1.sift')

figure()
gray()
plot_features(im1, l1, circle=True)
show()

（3）绘制出的特征点

7.SIFT特征匹配

（1）匹配代码如下

from PIL import Image
from pylab import *
import sys
from PCV.localdescriptors import sift


if len(sys.argv) >= 3:
  im1f, im2f = sys.argv[1], sys.argv[2]
else:
#  im1f = '../data/sf_view1.jpg'
#  im2f = '../data/sf_view2.jpg'
  im1f = 'D:\JMU\computer_vision\SIFTcode\SIFT\hle\JMU6.jpg'
  im2f = 'D:\JMU\computer_vision\SIFTcode\SIFT\hle\JMU14.jpg'
#  im1f = '../data/climbing_1_small.jpg'
#  im2f = '../data/climbing_2_small.jpg'
im1 = array(Image.open(im1f).convert("L"))
im2 = array(Image.open(im2f).convert("L"))
sift.process_image(im1f, 'out_sift_1.txt')
l1, d1 = sift.read_features_from_file('out_sift_1.txt')
figure()
gray()
subplot(121)
sift.plot_features(im1, l1, circle=False)

sift.process_image(im2f, 'out_sift_2.txt')
l2, d2 = sift.read_features_from_file('out_sift_2.txt')
subplot(122)
sift.plot_features(im2, l2, circle=False)

#matches = sift.match(d1, d2)
matches = sift.match_twosided(d1, d2)
print '{} matches'.format(len(matches.nonzero()[0]))

figure()
gray()
sift.plot_matches(im1, im2, l1, l2, matches, show_below=True)
show()

（2）匹配结果如下

                  

 

8.Harris特征匹配

（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("D:\JMU\computer_vision\SIFTcode\SIFT\hle\JMU3.jpg").convert("L"))
im2 = array(Image.open("D:\JMU\computer_vision\SIFTcode\SIFT\hle\JMU111.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()

（2）匹配结果如下
                  

9.多张图片索引出匹配度最高的

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


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 = loadtxt(filename)
    return f[:, :4], f[:, 4:]  # 特征位置，描述子


def write_featrues_to_file(filename, locs, desc):
    """将特征位置和描述子保存到文件中"""
    savetxt(filename, hstack((locs, desc)))


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

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


def match(desc1, desc2):
    """对于第一幅图像中的每个描述子，选取其在第二幅图像中的匹配
    输入：desc1(第一幅图像中的描述子)，desc2(第二幅图像中的描述子)"""
    desc1 = array([d / linalg.norm(d) for d in desc1])
    desc2 = array([d / linalg.norm(d) for d in desc2])
    dist_ratio = 0.6
    desc1_size = desc1.shape
    matchscores = zeros((desc1_size[0], 1), 'int')
    desc2t = desc2.T  # 预先计算矩阵转置
    for i in range(desc1_size[0]):
        dotprods = dot(desc1[i, :], desc2t)  # 向量点乘
        dotprods = 0.9999 * dotprods
        # 反余弦和反排序，返回第二幅图像中特征的索引
        indx = argsort(arccos(dotprods))
        # 检查最近邻的角度是否小于dist_ratio乘以第二近邻的角度
        if arccos(dotprods)[indx[0]] < dist_ratio * arccos(dotprods)[indx[1]]:
            matchscores[i] = int(indx[0])
    return matchscores


def match_twosided(desc1, desc2):
    """双向对称版本的match()"""
    matches_12 = match(desc1, desc2)
    matches_21 = match(desc2, 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):
    """返回将两幅图像并排拼接成的一幅新图像"""
    # 选取具有最少行数的图像，然后填充足够的空行
    rows1 = im1.shape[0]
    rows2 = im2.shape[0]
    if rows1 < rows2:
        im1 = concatenate((im1, zeros((rows2 - rows1, im1.shape[1]))), axis=0)
    elif rows1 > rows2:
        im2 = concatenate((im2, zeros((rows1 - rows2, im2.shape[1]))), axis=0)
    return 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 = 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')


# 获取project2_data文件夹下的图片文件名(包括后缀名)
filelist = get_imlist('D:\JMU\computer_vision\SIFTcode\SIFT')

# 输入的图片
im1f = 'D:\JMU\computer_vision\SIFTcode\SIFT\JMU111.jpg'
im1 = array(Image.open(im1f))
process_image(im1f, 'out_sift_1.txt')
l1, d1 = read_features_from_file('out_sift_1.txt')

i = 0
num = [0] * 30  # 存放匹配值
for infile in filelist:  # 对文件夹下的每张图片进行如下操作
    im2 = array(Image.open(infile))
    process_image(infile, 'out_sift_2.txt')
    l2, d2 = read_features_from_file('out_sift_2.txt')
    matches = match_twosided(d1, d2)
    num[i] = len(matches.nonzero()[0])
    i = i + 1
    print '{} matches'.format(num[i - 1])  # 输出匹配值

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

结果截图

 

10.实验分析

（1）通过SIFT算法和Harris算法对特征点的提取，我们可以发现使用SIFT算法提取到的特征点更多一些。因此我们也可以从此处知道SIFT算法更适用于图像特征检测，且比harris灵活、匹配准确些

（2）通过实验可以知道，SIFT算法的匹配正确性要高于Harris算法，因为由上面的图可知harris对旋转角度识别性不大，基本上出现匹配点错误，SIFT出现的较少。由上图也可知，sift通过主方向梯度改变会达到很好的旋转不变性效果。

11.遇到的问题与解决方法

（1）运行检测兴趣点代码时遇到如下问题：

通过百度查找资料，发现是因为找不到.sift文件，是由于VLfeat的版本问题。之后将版本换成了vlfeat-0.9.20的就好了。下载安装vlfeat包，版本0.9.20  官方下载连接：http://www.vlfeat.org/download/（注意要下载vlfeat-0.9.20版本）

（2）进行SIFT特征匹配时，出现如下问题

由此可知是因为输入的数组尺寸大小不一样导致的错误。解决方法：将两张要匹配的图片的大小及像素设置成一样的（或者是一同一宽度比列转换高度）。代码如下：

# -*- coding: utf-8 -*-
from PIL import Image
def produceImage(file_in, width, height, file_out):
    image = Image.open(file_in)
    resized_image = image.resize((width, height), Image.ANTIALIAS)
    resized_image.save(file_out)
 
if __name__ == '__main__':
    file_in = 'D:\JMU\computer_vision\SIFTcode\SIFT\hle\JMU111.jpg'
    width = 240
    height = 180
    file_out = 'D:\JMU\computer_vision\SIFTcode\SIFT\hle\JMU111new.jpg'
    produceImage(file_in, width, height, file_out)

通过上面代码，可以将图片转换成宽240，高180的图片。从而方便了图片的匹配