2021-03-28

计算机视觉 局部图像描述子

一、Harris角点检测

1.基本概念

确定特征点的位置，方向和尺度等仿射变换参数的过程。目前流行的特征点检测算法：Laplacian检测算法，高斯差分DOG算法，以及根据梯度协方差矩阵检测图中的角点位置的*等检测算法，其中基于DOG的特征点检测方法——尺度不变特征变换（sift）是目前性能最好的特征点检测算法。SIFT流程：对输入图像进行DOG滤波，搜索滤波后的图像中的所有的极大值和极小值，这些极值对应的像素坐标即为特征点坐标，当特征区域的大小与DOG滤波宽度大致相当时即呈现极值。

2.描述

SIFT提取特征点周围的一块区域，并利用特征向量来描述该区域作为特征描述符。SIFT描述符使用图像梯度，计算局部图像梯度方向的直方图，并在特征点周围创建4*4的直方图方格，每个柱状图包含八个梯度方向，因此SIFT 包含128维特征向量。
我们把图像域中点x上的对称半正定矩阵MI=MI(x)定义为：

其中∇I为包含导数Ix和Iy的图像梯度。由于该定义，MI的秩为1，特征值为λ1=|∇I|2和λ2=0。现在对于图像的每个像素，我们可以计算出该矩阵。
选择权重矩阵W（通常为高斯滤波器Gσ），我们可以得到卷积：

该卷积的目的是得到MI在周围像素上的局部平均。计算出的矩阵MI-有称为Harris矩阵。W的宽度决定了在像素x周围的感兴趣区域。像这样在区域附近对矩阵MI-取平均的原因是，特征值会依赖于局部图像特性而变化。如果图像梯度在该区域变化，那么MI-的第二个特征值将不再为0.如果图像的梯度没有变换，MI-的特征值也不会变化。
取决于该区域∇I的值，Harris矩阵MI-的特征值有三种情况：
1）如果λ1和λ2都是很大的正数，则该x点为角点；
2）如果λ1很大，λ2≈0，则该区域内存在一个边，该区域内的平均MI-的特征值不会变化太大；
3）如果λ1≈λ2≈0，则该区域为空。

为了去除加权常数κ，我们通常使用商数

作为指示器。

3.代码

from pylab import *
from PIL import Image
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 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')
    return

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
        # 反余弦和反排序，返回第二幅图像中特征的索引
        index = argsort(arccos(dotprods))

        # 检查最近邻的角度是否小于dist_ratio乘以第二近邻的角度
        if arccos(dotprods)[index[0]] < dist_ratio * arccos(dotprods)[index[1]]:
            matchscores[i] = int(index[0])

    return matchscores

def match_twosided(desc1,decs2):
    """双向对称版本的match"""

    matches_12 = match(desc1, decs2)
    matches_21 = match(decs2, decs2)

    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):
    """返回将两幅图像并排拼接成的一幅新图像"""

    # 选取具有最少行数的图像，然后填充足够的空行
    row1 = im1.shape[0]
    row2 = im2.shape[0]

    if row1 < row2:
        im1 = concatenate((im1,zeros((row2-row1,im1.shape[1]))), axis=0)
    elif row1 > row2:
        im2 = concatenate((im2,zeros((row1-row2,im2.shape[1]))), axis=0)

    # 如果这些情况都没有，那么他们的行数相同，不需要进行填充

    return concatenate((im1,im2), axis=1)

def plot_matches(im1: object, im2: object, locs1: object, locs2: object, matchscores: object, show_below: object = True) -> object:
    """显示一幅带有连接匹配之间连线的图片
        输入：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')

if __name__ == '__main__':
    # imname = 'raccoon.jpg'
    # im1 = array(Image.open(imname).convert('L'))
    # process_image(imname, 'raccoon.sift')
    # l1, d1 = read_features_from_file('raccoon.sift')
    #
    # figure()
    # gray()
    # plot_features(im1, l1, circle=True)
    # show()

    im1f = r'C:\Users\16089\Desktop\计算机视觉\cyk1.jpg'
    im2f = r'C:\Users\16089\Desktop\计算机视觉\cyk2.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()
    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()

4.结果

原图1：

原图2：

实验结果：

4.结论：

通过实验发现Harris角点的匹配数很多，SIFT匹配点数量少一些。
SIFT匹配点因为有通过当两幅图像的SIFT特征向量生成后，下一步我们采用关键点特征向量的欧式距离来作为两幅图像中关键点的相似性判定度量。取图像1中的某个关键点，并找出其与图像2中欧式距离最近的前两个关键点，在这两个关键点中，如果最近的距离除以次近的距离少于某个比例阈值，则接受这一对匹配点。降低这个比例阈值，SIFT匹配点数目会减少，但更加稳定。
非极大值抑制就是用一个模板，以该模板的中心点为准，计算该点的角点响应值R，如果该点的响应值R比周围的响应值都大，则周围的响应值变为0，如果该点的响应值R比周围的响应值都笑，则该点的响应值变为0

二、地理特征匹配

1.数据集

2.代码

# -*- coding: utf-8 -*-
import json
import os
import urllib
from pylab import *
from PIL import Image
from PCV.localdescriptors import sift
from PCV.tools import imtools
import pydot

if __name__ == '__main__':
    download_path = "D:\\image"
    path = "D:\\image"

    imlist = imtools.get_imlist(download_path)
    nbr_images = len(imlist)

    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

    # 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) + '.jpg'
                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) + '.jpg'
                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_jpg('22.jpg')

3.实验结果

3.实验总结

实验放置了三组图像：嘉庚图书馆，延奎图书馆，陆大楼，皆为不同视角所拍摄。
实验测试了7张图片，其中分别配对成功，虽然图片的角度与距离不太一样，还有部分障碍物遮挡，但仍然能匹配成功，所以说明SIFT算法在地理建筑匹配中的确有很大的优势。sift对旋转角度识别性比较大，而harris对旋转角度识别性不大，基本上出现匹配点错误，虽然sift也有匹配错误点，但是相对来说不会随意匹配，从harris匹配线可看出，匹配线杂乱无章，但是sift算子位置几乎是平行对应。所以对于旋转角度过后sift算法会更适合图像匹配，不过对于像素点多而又没有旋转角度拍摄的可以用harris来匹配，视觉效果可能会比较清晰，如果想匹配效果达到很高，建议还是使用sift算法匹配。