PCV库和vlfeat实现SIFT特征提取

1.原理

        1.1基本流程

                Lowe将SIFT算法分解为如下四步：

                1.尺度空间极值检测：搜索所有尺度上的图像位置。通过高斯差分函数来识别潜在的对于尺度和旋转不变的关键点。

                2.关键点定位：在每个候选的位置上，通过一个拟合精细的模型来确定位置和尺度。关键点的选择依据于它们的稳定程度。

                3.关键点方向确定：基于图像局部的梯度方向，分配给每个关键点位置一个或多个方向。所有后面的对图像数据的操作都相对于关键点的方向、尺度和位置进行变换，从而保证了对于这些变换的不变性。

                4.关键点描述：在每个关键点周围的邻域内，在选定的尺度上测量图像局部的梯度。这些梯度作为关键点的描述符，它允许比较大的局部形状的变形或光照变化。

我们就沿着Lowe的步骤，对SIFT算法的实现过程进行介绍：

        1.2尺度空间极值检测

                在不同的尺度空间是不能使用相同的窗口检测极值点，对小的关键点使用小的窗口，对大的关键点使用大的窗口，为了达到上述目的，我们使用尺度空间滤波器。

高斯核是唯一可以产生多尺度空间的核函数。-《Scale-space theory: A basic tool for analysing structures at different scales》。

        一个图像的尺度空间L(x,y,σ)，定义为原始图像I(x,y)与一个可变尺度的2维高斯函数G(x,y,σ)卷积运算。σ是尺度空间因子，它决定了图像的模糊的程度。在大尺度下（σ值大）表现的是图像的概貌信息，在小尺度下（σ值小）表现的是图像的细节信息。

        在计算高斯函数的离散近似时，在大概3σ距离之外的像素都可以看作不起作用，这些像素的计算也就可以忽略。所以，在实际应用中，只计算(6σ+1)*(6σ+1)的高斯卷积核就可以保证相关像素影响。

        下面我们构建图像的高斯金字塔，它采用高斯函数对图像进行模糊以及降采样处理得到的，高斯金字塔构建过程中，首先将图像扩大一倍，在扩大的图像的基础之上构建高斯金字塔，然后对该尺寸下图像进行高斯模糊，几幅模糊之后的图像集合构成了一个Octave，然后对该Octave下选择一幅图像进行下采样，长和宽分别缩短一倍，图像面积变为原来四分之一。这幅图像就是下一个Octave的初始图像，在初始图像的基础上完成属于这个Octave的高斯模糊处理，以此类推完成整个算法所需要的所有八度构建，这样这个高斯金字塔就构建出来了，整个流程如下图所示：

        利用LoG(高斯拉普拉斯方法)，即图像的二阶导数，可以在不同的尺度下检测图像的关键点信息，从而确定图像的特征点。但LoG的计算量大，效率低。所以我们通过两个相邻高斯尺度空间的图像的相减，得到DoG(高斯差分)来近似LoG。

        为了计算DoG我们构建高斯差分金字塔，该金字塔是在上述的高斯金字塔的基础上构建而成的，建立过程是：在高斯金字塔中每个Octave中相邻两层相减就构成了高斯差分金字塔。如下图所示：

        高斯差分金字塔的第1组第1层是由高斯金字塔的第1组第2层减第1组第1层得到的。以此类推，逐组逐层生成每一个差分图像，所有差分图像构成差分金字塔。概括为DOG金字塔的第o组第l层图像是有高斯金字塔的第o组第l+1层减第o组第l层得到的。后续Sift特征点的提取都是在DOG金字塔上进行的

        在 DoG 搞定之后，就可以在不同的尺度空间中搜索局部最大值了。对于图像中的一个像素点而言，它需要与自己周围的 8 邻域，以及尺度空间中上下两层中的相邻的 18（2x9）个点相比。如果是局部最大值，它就可能是一个关键点。基本上来说关键点是图像在相应尺度空间中的最好代表。如下图所示：

        搜索过程从每组的第二层开始，以第二层为当前层，对第二层的DoG图像中的每个点取一个3×3的立方体，立方体上下层为第一层与第三层。这样，搜索得到的极值点既有位置坐标（DoG的图像坐标），又有空间尺度坐标（层坐标）。当第二层搜索完成后，再以第三层作为当前层，其过程与第二层的搜索类似。当S=3时，每组里面要搜索3层，所以在DOG中就有S+2层，在初使构建的金字塔中每组有S+3层。 

        1.3关键点定位

        由于DoG对噪声和边缘比较敏感，因此在上面高斯差分金字塔中检测到的局部极值点需经过进一步的检验才能精确定位为特征点。

        使用尺度空间的泰勒级数展开来获得极值的准确位置， 如果极值点的灰度值小于阈值（一般为0.03或0.04）就会被忽略掉。 在 OpenCV 中这种阈值被称为 contrastThreshold。

        DoG 算法对边界非常敏感， 所以我们必须要把边界去除。 Harris 算法除了可以用于角点检测之外还可以用于检测边界。从 Harris 角点检测的算法中，当一个特征值远远大于另外一个特征值时检测到的是边界。那在DoG算法中欠佳的关键点在平行边缘的方向有较大的主曲率，而在垂直于边缘的方向有较小的曲率，两者的比值如果高于某个阈值（在OpenCV中叫做边界阈值），就认为该关键点为边界，将被忽略，一般将该阈值设置为10。

        将低对比度和边界的关键点去除，得到的就是我们感兴趣的关键点。

        1.4关键点方向确定

        经过上述两个步骤，图像的关键点就完全找到了，这些关键点具有尺度不变性。为了实现旋转不变性，还需要为每个关键点分配一个方向角度，也就是根据检测到的关键点所在高斯尺度图像的邻域结构中求得一个方向基准。

        对于任一关键点，我们采集其所在高斯金字塔图像以r为半径的区域内所有像素的梯度特征（幅值和幅角），半径r为：

 r=3×1.5σ

其中σ是关键点所在octave的图像的尺度，可以得到对应的尺度图像。

        梯度的幅值和方向的计算公式为：

        邻域像素梯度的计算结果如下图所示：

        完成关键点梯度计算后，使用直方图统计关键点邻域内像素的梯度幅值和方向。具体做法是，将360°分为36柱，每10°为一柱，然后在以r为半径的区域内，将梯度方向在某一个柱内的像素找出来，然后将他们的幅值相加在一起作为柱的高度。因为在r为半径的区域内像素的梯度幅值对中心像素的贡献是不同的，因此还需要对幅值进行加权处理，采用高斯加权，方差为1.5σ。如下图所示，为简化图中只画了8个方向的直方图。 

        每个特征点必须分配一个主方向，还需要一个或多个辅方向，增加辅方向的目的是为了增强图像匹配的鲁棒性。辅方向的定义是，当一个柱体的高度大于主方向柱体高度的80%时，则该柱体所代表的的方向就是给特征点的辅方向。

        直方图的峰值，即最高的柱代表的方向是特征点邻域范围内图像梯度的主方向，但该柱体代表的角度是一个范围，所以我们还要对离散的直方图进行插值拟合，以得到更精确的方向角度值。利用抛物线对离散的直方图进行拟合，如下图所示：

        获得图像关键点主方向后，每个关键点有三个信息(x,y,σ,θ)：位置、尺度、方向。由此我们可以确定一个SIFT特征区域。通常使用一个带箭头的圆或直接使用箭头表示SIFT区域的三个值：中心表示特征点位置，半径表示关键点尺度，箭头表示方向。如下图所示：

         1.5关键点描述

        通过以上步骤，每个关键点就被分配了位置，尺度和方向信息。接下来我们为每个关键点建立一个描述符，该描述符既具有可区分性，又具有对某些变量的不变性，如光照，视角等。而且描述符不仅仅包含关键点，也包括关键点周围对其有贡献的的像素点。主要思路就是通过将关键点周围图像区域分块，计算块内的梯度直方图，生成具有特征向量，对图像信息进行抽象。

        描述符与特征点所在的尺度有关，所以我们在关键点所在的高斯尺度图像上生成对应的描述符。以特征点为中心，将其附近邻域划分为d∗d个子区域（一般取d=4)，每个子区域都是一个正方形，边长为3σ，考虑到实际计算时，需进行三次线性插值，所以特征点邻域的为3σ(d+1)∗3σ(d+1)的范围，如下图所示：

        为了保证特征点的旋转不变性，以特征点为中心，将坐标轴旋转为关键点的主方向，如下图所示： 

        计算子区域内的像素的梯度，并按照σ=0.5d进行高斯加权，然后插值计算得到每个种子点的八个方向的梯度，插值方法如下图所示：

        每个种子点的梯度都是由覆盖其的4个子区域插值而得的。如图中的红色点，落在第0行和第1行之间，对这两行都有贡献。对第0行第3列种子点的贡献因子为dr，对第1行第3列的贡献因子为1-dr，同理，对邻近两列的贡献因子为dc和1-dc，对邻近两个方向的贡献因子为do和1-do。则最终累加在每个方向上的梯度大小为：

         其中k，m，n为0或为1。 如上统计4∗4∗8=1284∗4∗8=128个梯度信息即为该关键点的特征向量，按照特征点的对每个关键点的特征向量进行排序，就得到了SIFT特征描述向量。

2.SIFT特征提取实现 

        2.1安装PCV库

                2.1.1下载PCV包

                下载链接：https://codeload.github.com/Li-Shu14/PCV/zip/master；不用解压，copy到项目文件夹下；

                2.1.2安装PCV

                项目文件夹下运行cmd，conda activate进入指定虚拟环境，执行以下命令安装：

        pip install PCV-master.zip

​        2.2安装vlfeat

                1.下载地址：VLFeat - Home

                2.下载0.9.20版本，即vlfeat-0.9.20-bin.tar.gz

                 3.解压后进入bin文件夹，根据系统选择对应版本，很重要！！！将对应版本文件夹复制到项目文件夹下，这边将其重新命名为win64vlfeat;

                4.找到PCV包下localdescriptors文件夹下的sift.py文件，找到下面这段代码

def process_image(imagename,resultname,params="--edge-thresh 10 --peak-thresh 5"):
    """ Process an image and save the results in a file. """

    if imagename[-3:] != 'pgm':
        # create a pgm file
        im = Image.open(imagename).convert('L')
        im.save('tmp.pgm')
        imagename = 'tmp.pgm'

    cmmd = str(r"H:\2023.4sift\win64vlfeat\sift.exe "+imagename+" --output="+resultname+
                " "+params)
    os.system(cmmd)
    print('processed', imagename, 'to', resultname)

                将其中cmmd后的sift.exe文件路径改为上一步路径，注意路径前加r，sift.exe后有空格

        2.3描述子提取

from PIL import Image
#from numpy import *
from pylab import *
import os
import numpy as np
def process_image(imagename,resultname,params="--edge-thresh 10 --peak-thresh 5"):
    """ Process an image and save the results in a file. """
 
    if imagename[-3:] != 'pgm':
        # create a pgm file
        im = Image.open(imagename).convert('L')  #.convert('L') 将RGB图像转为灰度模式，灰度值范围[0,255]
        im.save('tmp.pgm')                       #将灰度值图像信息保存在.pgm文件中
        imagename = 'tmp.pgm'
   
    cmmd = str(r"H:/2023.4sift/win64vlfeat/sift.exe "+imagename+" --output="+resultname+
                " "+params)
    os.system(cmmd)                              #执行sift可执行程序，生成resultname(test.sift)文件
    print('processed', imagename, 'to', resultname)
 

#从上面输出文件中，将特征读取到Numpy数组中的函数
def read_features_from_file(filename):
    """ Read feature properties and return in matrix form. """
    #print(filename)
    f = np.loadtxt(filename)
    return f[:,:4],f[:,4:] # feature locations, descriptors
 
#读取特征后，通过在图像上绘制出它们的位置，可以将其可视化 
def plot_features(im,locs,circle=True):
    """ Show image with features. input: im (image as array), 
        locs (row, col, scale, orientation of each feature). """
 
    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')
 
 
if __name__ == '__main__':
    imname = r'H:\2023.4sift\1.jpg'       #待处理图像路径
    im=Image.open(imname)
    process_image(imname,'out_sift_1.txt')
    l1,d1 = read_features_from_file('out_sift_1.txt')  
    figure()
    gray()
    plot_features(im,l1,circle = True)
    title('sift-features')
    show()

原图：

运行结果：

 

得到的out_sift_1.txt文件:

 前四个数值分别为坐标、尺度和方向角度，后面跟对应描述符的128维向量。

         2.4配对描述子

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))
        
    plt.figure(figsize=(20, 10))
    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')

 结合2.3完整代码如下：

from PIL import Image
# from numpy import *
from pylab import *
import os
import numpy as np


def process_image(imagename, resultname, params="--edge-thresh 10 --peak-thresh 5"):
    """ Process an image and save the results in a file. """

    if imagename[-3:] != 'pgm':
        # create a pgm file
        im = Image.open(imagename).convert('L')  # .convert('L') 将RGB图像转为灰度模式，灰度值范围[0,255]
        im.save('tmp.pgm')  # 将灰度值图像信息保存在.pgm文件中
        imagename = 'tmp.pgm'

    cmmd = str(r"H:/2023.4sift/win64vlfeat/sift.exe " + imagename + " --output=" + resultname +
               " " + params)
    os.system(cmmd)  # 执行sift可执行程序，生成resultname(test.sift)文件
    print('processed', imagename, 'to', resultname)


def read_features_from_file(filename):
    """ Read feature properties and return in matrix form. """
    # print(filename)
    f = np.loadtxt(filename)
    return f[:, :4], f[:, 4:]  # feature locations, descriptors


def plot_features(im, locs, circle=True):
    """ Show image with features. input: im (image as array),
        locs (row, col, scale, orientation of each feature). """

    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))

    plt.figure(figsize=(20, 10))
    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__':
    # imname = (r'H:\2023.4sift\2.jpg')  # 待处理图像路径
    # im = Image.open(imname)
    # process_image(imname, 'out_sift_2.txt')
    # l1, d1 = read_features_from_file('out_sift_2.txt')  
    # figure()
    # gray()
    # plot_features(im, l1, circle=True)
    # title('sift-features')
    # show()
    im1f = r'H:\2023.4sift\1.jpg'
    im2f = r'H:\2023.4sift\2.jpg'
    im1 = array(Image.open(im1f))
    im2 = array(Image.open(im2f))

    # process_image(im1f, 'out_sift_1.txt')
    l1, d1 = read_features_from_file('out_sift_1.txt')
    figure()
    gray()
    plt.figure(figsize=(20, 10))
    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()

运行结果：

 参考链接：计算机视觉3 SIFT特征提取与全景图像拼接_print 'processed', imagename, 'to', resultname_吃花椒的恩酱的博客-CSDN博客

SIFT算法原理：SIFT算法详细介绍