[游戏开发]图像局部描述符【计算机视觉】

一、Harris角点检测算法

1.什么是角点

下面有两幅不同视角的图像，通过找出对应的角点进行匹配

我们可以直观的概括下角点所具有的特征：

轮廓之间的交点；

对于同一场景，即使视角发生变化，通常具备稳定性质的特征；

该点附近区域的像素点无论在梯度方向上还是其梯度幅值上有着较大变化；

2. 角点检测算法基本思想是什么？

移动窗口W，位移为（u,v），比较移动前后的像素变化
其中窗口函数可以是平坦的，也可以是高斯的如下图：

首先，将图像窗口平移[u,v]产生灰度变化的自相关函数如下：

经过一系列E(u，v)表达式的演化， E(u,v)表达式可以更新为：

其中矩阵M是2x2矩阵，可由图像的导数求得:

3.角点响应函数R

其中k是常量，一般取值为0.04~0.06，这个参数仅仅是这个函数的一个系数，它的存在只是调节函数的形状而已。

4. 用python实现Harris角点检测

（1）Harris角点检测器的响应函数

Harris角点检测器的响应函数会返回像素值为 Harris 响应函数值的一幅图像。

from PIL import Image
from numpy import *
from pylab import *
from scipy.ndimage import filters

def compute_harris_response(im,sigma=3): #在灰度图像中计算每个像素的Harris的焦点响应函数
    # 计算导数
    imx = zeros(im.shape)
    # 高斯倒数
    filters.gaussian_filter(im,(sigma,sigma),(0,1),imx)
    imy = zeros(im.shape)
    filters.gaussian_filter(im,(sigma,sigma),(1,0),imy)

    #矩阵分量
    Wxx = filters.gaussian_filter(imx*imx,sigma)
    Wxy = filters.gaussian_filter(imx*imy,sigma)
    Wyy = filters.gaussian_filter(imy*imy,sigma)

    #特征值
    Wdet = Wxx*Wyy - Wxy**2
    Wtr = Wxx + Wyy

    return Wdet/Wtr

(2)返回Harris角点函数

def get_harris_points(harrisim,min_dist=10,threshold=0.5):
    #挑选高于阈值
    corner_threshold = threshold
    harrisim_t = (harrisim > corner_threshold) * 1

    #得到坐标
    coords = array(harrisim_t.nonzero()).T
    #响应值
    candidate_values = [harrisim[c[0],c[1]] for c in coords]

    #按照响应值排序
    index = argsort(candidate_values)[::-1]

    #将可行点的位置存放于数组
    allowed_locations = zeros(harrisim.shape)
    allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1

    #选择最佳的角点
    filtered_coords = []
    for i in index:
        if allowed_locations[coords[i,0],coords[i,1]] == 1:
            filtered_coords.append(coords[i])
            allowed_locations[(coords[i,0]-min_dist):(coords[i,0]+min_dist),(coords[i,1]-min_dist):(coords[i,1]+min_dist)] = 0
    return filtered_coords

（3）可视化

def plot_harris_points(image,filtered_coords):
    #可视化
    figure(dpi = 120)
    gray()
    imshow(image)
    plot([p[1] for p in filtered_coords],
         [p[0] for p in filtered_coords],'*')
    axis('off')
    show()

（4）Harris角点检测实例

from pylab import *
from PIL import Image

#读入图像
im = array(Image.open('data/test/no/jmu1.jpg').convert('L'))
#检测harris响应函数
harrisim = compute_harris_response(im)
#Harris响应函数
harrisiml = 255 - harrisim

figure(dpi = 150)
gray()
#画出响应图
subplot(221)
imshow(im)
axis('off')
axis('equal')

threshold = [0.01,0.05,0.1]
for i, thres in enumerate(threshold):
    filtered_coords = get_harris_points(harrisim,50,thres)
    subplot(2,2,i+2)
    title(str(threshold[i]))
    imshow(im)
    plot([p[1] for p in filtered_coords],[p[0] for p in filtered_coords], '*')
    axis('off')

show()

（5）结果展示

5. 用opencv实现Harris角点检测

API 在opencv中有提供实现 Harris 角点检测的函数 cv2.cornerHarris，我们直接调用的就可以，非常方便。

函数原型：cv2.cornerHarris(src, blockSize, ksize, k[, dst[, borderType]])

import cv2
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
np.seterr (divide='ignore',invalid='ignore')

filename = 'data/test/no/123.jpg'
img = cv2.imread(filename)
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

gray = np.float32(gray)
dst = cv2.cornerHarris(gray,2,3,0.04)

img[dst > 0.01 * dst.max()] = [0,255,0]

plt.figure(dpi = 120)
plt.imshow(img[: ,: ,[2,1,0]])
plt.show()

结果展示

6. Harris角点检测算法优缺点：

（1）优点

1.旋转不变性，椭圆转过一定角度但是其形状保持不变（特征值保持不变）。
2.对于图像灰度的仿射变化具有部分的不变性，由于仅仅使用了图像的一介导数，对于图像灰度平移变化不变；
3.对于图像灰度尺度变化不变。

（2）缺点：

1.它对尺度很敏感，不具备几何尺度不变性。
2.提取的角点是像素级的

二、SIFT（尺度不变特征变换）

1.SIFT概述

SIFT，即尺度不变特征变换（Scale-invariant feature transform，SIFT），是用于图像处理领域的一种描述。这种描述具有尺度不变性，它对物体的尺度变化，刚体变换，光照强度和遮挡都具有较好的稳定性，可在图像中检测出关键点，是一种局部特征描述子。SIFT 算法被认为是图像匹配效果好的方法之 一算法实现特征匹配主要有三个流程：

①特征点提取；
②特征点主方向确定；
③特征点描述；
④特征点匹配；

其中特征点提取主要包括生成高斯差分（DifferenceofGaussian，DOG）尺度空间、寻找局部极值点、特征点筛选、确定特征点方向；特征点匹配主要包括根据描述子相似性进行匹配、匹配对比值提纯、RANSAC方法剔除离群匹配对。

2.SIFT特征提取及匹配

（1）特征点提取（关键点检测）

两种图像在匹配的时候可能因为拍摄的距离、拍摄的角度问题，会导致在特征点提取的时候差异很大，所以我们希望SIFT的特征点可以具有尺度不变性和方向不变性。

（2） 尺度空间的构建

图像的尺度空间是这幅图像在不同解析度下的表示。一幅图像可以产生几组（octave）图像，一组图像包括几层图像。构造尺度空间传统的方法即构造一个高斯金字塔，原始图像作为最底层，然后对图像进行高斯模糊再降采样（2倍）作为下一层图像（即尺度越大，图像越模糊），循环迭代下去。

对图像进行尺度变换，以满足特征点的尺度不变性，保留图像轮廓和细节。

(3)DOG算子（Difference of Gaussian）

DoG（Difference of Gaussian）函数：

该函数在计算上只需相邻高斯平滑后图像相减，因此简化了计算。

（4）DOG局部极值检测

特征点是由DOG空间的局部极值点组成的。为了寻找DoG函数的极值点， 每一个像素点要和它所有的相邻点比较，看其是否比它的图像域和尺度域的相邻点大或者小。

中间的检测点和它同尺度的8个相邻点和上下相邻尺度对应的9×2个 点共26个点比较，以确保在尺度空间和二维图像空间都检测到极值点。

有些极值点的位置是在图像的边缘位置的，因为图像的边缘点很难定位，同时也容易受到噪声的干扰，我们把这些点看做是不稳定的极值点，需要进行去除。边缘梯度的方向上主曲率值比较大，而沿着边缘方向则主曲率值较小。候选特征点的DoG函数D(x)的主曲率与2×2Hessian矩阵H的特征值成正比，Dxx表示DOG金字塔中某一尺度的图像x方向求导两次：

在在边缘梯度的方向上主曲率值比较大，而沿着边缘方向则主曲率值较小。设α=λmax=Dxx为H的最大特征值,β=λmin=Dyy为H的最小特征值，则有：

Tr(H) 为矩阵H的迹,Det(H)为矩阵H的行列式.当两个特征值相等时其值最小,因此为了检测主曲率是否在某个阈值Tr下，只需检测该比值与阈值T的大小关系，过滤不稳定的边缘响应点。

所以，特征点提取可以概括为以下几个步骤：

1.构建高斯尺度空间，产生不同尺度的高斯模糊图像。
2.进行降采用，得到一系列尺寸不断缩小的图像。
3.DOG空间极值检测，去除部分边缘响应点。

3.特征点描述

通过以上步骤，对于每一个关键点，拥有三个信息：位置、尺度以及方向。接下来就是为每个关键点建立一个描述符，使其不随各种变化而改变，比如光照变化、视角变化等等。并且描述符应该有较高的独特性，以便于提高特征点正确匹配的概率。

4.如何实现特征匹配?

在不同尺度空间上查找特征点（关键点）的问题)
1、提取关键点;
2、对关键点附加详细的信息(局部特征)，即描述符;
3、通过特征点(附带上特征向量的关键点)的两两比较找出相互匹配的若干对特征点，建立景物间的对应关系

5.用python检测兴趣点

（1）将图像转换成.pgm格式文件

def process_image(imagename,resultname,params="--edge-thresh 10 --peak-thresh 5"):



    if imagename[-3:] != 'pgm':
        im = Image.open(imagename).convert('L')
        im.save('tmp.pgm')
        imagename = 'tmp.pgm'
    cmmd = str(r'D:\SVMCV\c02\sift.exe '+imagename+" --output="+resultname+" "+params)
    os.system(cmmd)
    print('processed', imagename,'to',resultname)

（2）读取图像特征值属性值

读取特征值属性值，然后将其以矩阵形式返回

def read_features_from_file(filename):
    f = loadtxt(filename)
    return f[:,:4],f[:,4:]

（3） 可视化

def plot_features(im,locs,circle=False):
    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

（4）检测兴趣点实例并与Harris角点检测对比

imname = 'data/test/no/123.jpg'

im = array(Image.open(imname).convert('L'))
#得到SIFT特征
process_image(imname,'122.sift')
#取出坐标以及特征描述
l1,d1 = read_features_from_file('122.sift')
figure(dpi=180)
gray()
subplot(221)
plt.axis('off')
imshow(Image.open(imname))
title('src')
subplot(222)
#可视化特征点
plot_features(im,l1,circle=False)
title('SIFT')
subplot(223)
plot_features(im,l1,circle=True)
title('scale')
#检测harris角点，对比
harrisim = compute_harris_response(im)
subplot(224)
filtered_coords = get_harris_points(harrisim,60,0.1)
imshow(im)
plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '*')
axis('off')
title('Harris')
show()

（5）结果展示

（6）结果分析

sift和Harris角点的两种算子选择了不同的坐标
SIFT算法的实质是在不同的尺度空间上查找关键点(特征点)，并计算出关键点的方向。SIFT所查找到的关键点是一些十分突出，不会因光照，仿射变换和噪音等因素而变化的点，如角点、边缘点、暗区的亮点及亮区的暗点等。

6.用SIFT实现特征匹配

对于将一幅图像中的特征匹配到另一幅图像的特征，一种稳健的准则（同样是由Lowe提出的）是使用者两个特征距离和两个最匹配特征距离的比率。相比于图像中的其他特征，该准则保证能够找到足够相似的唯一特征。使用该方法可以使错误的匹配数降低。

代码展示

im1f = 'data/test/no/123.jpg'
im2f = 'data/test/no/321.jpg'

#sift特征提取
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)



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


    # 选取具有最少行数的图像，然后填充足够的空行
    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, 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[int(matchscores[i, 0]), 0] + cols1], [locs1[i, 1], locs2[matchscores[i, 0], 1]], 'c')
    axis('off')

#特征匹配
matches = match_twosided(d1,d2)
print ('{} matches'.format(len(matches.nonzero()[0])))
figure(dpi = 180)
gray()
plot_matches(im1,im2,l1,l2,matches,show_below=True)
show()

结果展示

三，匹配地理标记图像

代码

# -*- codeing =utf-8 -*-
# @Time : 2021/3/23 10:45
# @Author : ArLin
# @File : demo8.py
# @Software: PyCharm
# -*- 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 = "D:\imagepi"  # set this to the path where you downloaded the panoramio images
path = "D:\imagepi"  # 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((200, 200))
            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('llJMU.png')

结果展示