计算机视觉：特征提取与匹配

目录1. 特征提取和匹配1.1 背景知识1.2 特征匹配基本流程1.3局部特征描述子2. Harris角点检测2.1 角点（corner points）2.2HARRIS角点检测基本思想2.3HARRIS检测：数学表达2.4 角点响应函数2.5 编程实现2.5.1 角点检测代码实现2.5.2 角点检测数据测试3.SIFT特征匹配算法3.1 SIFT算法实现步骤3.2关键点检测的相关概念3.2.1

DreamWendy

23775人浏览 · 2022-03-31 01:14:25

DreamWendy · 2022-03-31 01:14:25 发布

2.1 角点（corner points）

3.3.1 DoG（Difference of Gaussian）函数

1. 特征提取和匹配

1.1 背景知识

特征提取和匹配是许多计算机视觉应用中的一个重要任务，广泛运用在运动结构、图像检索、目标检测等领域。特征可能是图像中的特定结构，例如点，边缘或对象。特征也可能是应用于图像的一般邻域操作或特征检测的结果。特征点匹配，是图像拼接、三维重建、相机标定等应用的关键步骤。

关键点也称兴趣点，是纹理中表达的点。关键点往往是物体边界方向突然改变的点或两个或多个边缘段之间的交点。它在图像空间中具有明确的位置或很好地定位。即使图像域的局部或全局存在如光照和亮度变化等的扰动，关键点仍然是稳定，可以被重复可靠地计算出。除此之外它能够提供有效的检测。

1.2 特征匹配基本流程

根据特定准则，提取图像中的特征点；
提取特征点周围的图像块，构造特征描述符；
通过特征描述符对比，实现特征匹配。

1.3 局部特征描述子

特征点必须具有不变性，比如：

几何不变性: 位移、旋转、尺度，......
光度不变性: 光照、曝光，......

特征分类

角点。Harris算子，SUSAN算子, FAST算子。
梯度特征点。 SIFT、SURF、GLOH、ASIFT、PSIFT算子等。  边缘特征（线型）。Canny算子, Marr算子。
纹理特征。灰度共生矩阵，小波Gabor算子。
快速特征匹配的基本策略：以图像中稳定角点的领域训练特征分类器，为输入图像的特征分类。

2. Harris角点检测

Harris 角点检测算法（也称 Harris & Stephens 角点检测器）是一个极为简单的角点检测算法。该算法的主要思想是，如果像素周围显示存在多于一个方向的边，我们认为该点为兴趣点。该点就称为角点。

2.1 角点（corner points）

局部窗口沿各方向移动，均产生明显变化的点；
图像局部曲率突变的点。

好的角点检测算法应该具有的特点：

检测出图像中“真实的”角点
准确的定位性能
很高的稳定性
具有对噪声的鲁棒性
具有较高的计算效率

2.2 HARRIS角点检测基本思想

算法基本思想是使用一个固定窗口在图像上进行任意方向上的滑动，比较滑动前与滑动后两种情况，窗口中的像素灰度变化程度，如果存在任意方向上的滑动，都有着较大灰度变化，那么我们可以认为该窗口中存在角点。

固定窗口在图像上进行任意方向上的滑动，有以下几种情况：

平坦区域：任意方向移动，无灰度变化。

平坦区域

边缘：沿着边缘方向移动，无灰度变化。

边缘上移动

角点：窗口向任意方向的移动都导致图像灰度的明显变化。

角点处移动

2.3 HARRIS检测：数学表达

将图像窗口平移 $[u,v]$ 产生灰度变化 $E(u,v)$ 的数学表达式如下：

$E(u,v)=\sum_{x,y}^{}w(x,y)[I(x+u,y+v)-I(x,y)]^2$

其中， $w(x,y)$ 为窗口函数，可用1表示在窗口内部，用0表示在外面，或者用高斯拉普拉斯算子。

$I(x,y)$ 为灰度图像， $I(x+u,y+v)$ 为平移后的图像灰度。

延伸知识点：

一元函数泰勒展开：

二元函数泰勒展开：

为了对 $I(x+u,y+v),E(u,v)$ 求解方便，我们利用二元泰勒近似可将表达式简化，将 $I(x+u,y+v)$ 函数在 $(x,y)$ 处泰勒展开，得：

$I(x+u,y+v)=I(x,y)+I_xu+I_yv+O(x^2,y^2)$

$E(u,v)=\sum_{x,y}^{}w(x,y)[I_xu+I_yv+O(x^2,y^2)]^2 \approx \sum_{x,y}^{}w(x,y)[I_xu+I_yv]^2$

又我们可以将 $[I_xu+I_yv]^2$ 表示为：

$[I_xu+I_yv]^2=[u,v]\begin{bmatrix} I_x^2 &I_xI_y \\ I_xI_y& I_y^2 \end{bmatrix}\begin{bmatrix} u\\ v \end{bmatrix}$

于是对于局部微小的移动量 $[u,v]$ ，可以近似得到下面的表达:

$E(u,v)\cong [u,v]M\begin{bmatrix} u\\ v \end{bmatrix}$

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

$M=\sum_{x,y}^{}w(x,y)\begin{bmatrix} I_x^2 &I_xI_y \\ I_xI_y& I_y^2 \end{bmatrix}$

窗口移动导致的图像变化量：实对称矩阵M的特征值分析

对于， $E(u,v)\cong [u,v]M\begin{bmatrix} u\\ v \end{bmatrix}$ ，记 $M$ 的特征值为 $\lambda_1, \lambda_2$ .

通过M的两个特征值的大小对图像点进行分类:

平台区域：如果 $\lambda_1$ 和 $\lambda_2$ 都很小，图像窗口在所有方向上移动都无明显灰度变化;
边缘：如果 $\lambda_1>> \lambda_2$ 或 $\lambda_1<< \lambda_2$ ，图像窗口在某一方向上移动有明显灰度变化。
角点：如果 $\lambda_1$ 和 $\lambda_2$ 都较大且数值相当 $\lambda_1$ ~ $\lambda_2$ ，图像窗口在所有方向上移动都产生明显灰度变化;

2.4 角点响应函数

在不需要实际计算特征值的情况下，为了把重要的情况和其他情况分开，Harris 和 Stephens 引入了指示函数。

定义：角点响应函数R

$R=det M-k(trace M)^2$

$det M=\lambda_1 \lambda _2$

$trace M=\lambda_1 +\lambda _2$

$(k - empirical constant, k = 0.04\sim 0.06)$

R 只与M的特征值有关

角点：R 为大数值正数
边缘：R为大数值负数
平坦区：R为小数值

角点计算流程

对角点响应函数R进行阈值处理：R>threshold

提取R的局部最大值

为了消除参数k的影响，也可采用商来计算响应：

$R=\frac{det M}{(trace M)^2}$

$det M=\lambda_1 \lambda _2$

$trace M=\lambda_1 +\lambda _2$

2.5 编程实现

Harris角点检测步骤

读入图像转灰度图

计算角点响应函数

基于响应值选择角点

在原点中画出检测的角点

2.5.1 角点检测代码实现

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

# 在灰度图像中，对每个像素计算Harris角点检测器响应函数
def compute_harris_response(im, sigma=3):

    # 计算导数: 先初始化矩阵
    imx = zeros(im.shape)
    imy = zeros(im.shape)
    # 高斯导数: 3*3的高斯滤波
    filters.gaussian_filter(im, (sigma, sigma), (0, 1), imx)
    filters.gaussian_filter(im, (sigma, sigma), (1, 0), imy)

    # 计算 Harris 矩阵的分量
    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

# 从一幅Harris响应图像中返回角点
def get_harris_points(harrisim, min_dist=10, threshold=0.1):
    # min_dist 为分割角点和图像边界的最少像素数目

    # 寻找高于阈值的候选角点
    corner_threshold = harrisim.max() * threshold
    harrisim_t = (harrisim > corner_threshold) * 1

    # 得到候选点的坐标
    coords = array(harrisim_t.nonzero()).T
    # 获取候选点的 Harris 响应值
    candidate_values = [harrisim[c[0], c[1]] for c in coords]
    # 对候选点按照 Harris 响应值进行排序
    index = argsort(candidate_values)
    # 将可行点的位置保存到数组中
    allowed_locations = zeros(harrisim.shape)
    allowed_locations[min_dist:-min_dist, min_dist:-min_dist] = 1
    # 按照 min_distance 原则，选择最佳 Harris 点
    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

# 绘制图像中检测到的角点
def plot_harris_points(image, filtered_coords):
    figure()
    gray()
    imshow(image)
    plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '*')
    axis('off')
    show()

2.5.2 角点检测数据测试

数据一：

数据二：

数据三：

测试代码示例：

if __name__ == '__main__':

    # 读入图像
    im = array(Image.open('Images/picture1.jpg').convert('L'))
    # 检测Harris角点
    harrisim = compute_harris_response(im)
    # Harris响应函数
    harrisim1 = 255 - harrisim

    figure()
    gray()
    # 画出Harris响应图
    subplot(221)
    title("harrisim")
    imshow(harrisim1)
    axis('off')
    axis('equal')

    # 设置多个阈值，比较实验结果
    threshold = [0.1, 0.2, 0.3]
    for i, thres in enumerate(threshold):
        filtered_coords = get_harris_points(harrisim, 6, thres)
        subplot(2, 2, i+2)
        title("threshold: "+str(threshold[i]))
        imshow(im)
        plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '*')
        axis('off')
    show()

运行结果：

数据一结果：

数据二结果：

数据三结果：

阈值增大，角点响应值R减小，降低角点检测的灵性，被检测角点数量减少；阈值减小，角点响应值R增大，增加角点检测的灵敏性，被检测角点的数量增加。

3. SIFT特征匹配算法

David Lowe 提出的 SIFT（Scale-Invariant Feature Transform，尺度不变特征变换）是过去十年中最成功的图像局部描述子之一。SIFT 特征后来得到精炼并详述，经受住了时间的考验。SIFT 特征包括兴趣点检测器和描述子。SIFT 描述子具有非常强的稳健性，这在很大程度上也是 SIFT 特征能够成功和流行的主要原因。自从 SIFT 特征的出现，许多其他本质上使用相同描述子的方法也相继出现。现在，SIFT 描述符经常和许多不同的兴趣点检测器相结合使用（有些情况下是区域检测器），有时甚至在整幅图像上密集地使用。SIFT 特征对于尺度、旋转和亮度都具有不变性，因此，它可以用于三维视角和噪声的可靠匹配。

SIFT算法可以解决的问题：

目标的旋转、缩放、平移（RST）
图像仿射/投影变换（视点viewpoint）
弱光照影响（illumination）
部分目标遮挡（occlusion）
杂物场景（clutter）
噪声

3.1 SIFT算法实现步骤

实质可以归为在不同尺度空间上查找特征点（关键点）的问题。

SIFT算法实现特征匹配主要有三个流程：

提取关键点；
对关键点附加详细的信息（局部特征），即描述符；
通过特征点（附带上特征向量的关键点）的两两比较找出相互匹配的若干对特征点，建立景物间的对应关系。

3.2 关键点检测的相关概念

3.2.1 关键点（特征点）

这些点是一些十分突出的点,不会因光照、尺度、旋转等因素的改变而消失，比如角点、边缘点、暗区域的亮点以及亮区域的暗点。假定两幅图像中有相同的景物，那么使用某种算法分别提取各自的特征点，这些点之间会有相互对应的匹配关系。根据归纳，我们可以看出SIFT特征点希望选出具有下述不变性的点：尺度、方向、位移、光照。

3.2.2 尺度空间

尺度空间理论最早于1962年提出，其主要思想是通过对原始图像进行尺度变换，获得图像多尺度下的空间表示。从而实现边缘、角点检测和不同分辨率上的特征提取，以满足特征点的尺度不变性。

尺度空间中各尺度图像的模糊程度逐渐变大，能够模拟人在距离目标由近到远时目标在视网膜上的形成过程。尺度越大图像越模糊。

尺度不变性

旋转不变性

根据文献《Scale-space theory: A basic tool for analysing structures at different scales》可知，高斯核可以产生多尺度空间的核，一个图像的尺度空间L(x, y, σ) ，定义为原始图像 I(x, y)与一个可变尺度的2维高斯函数G(x, y, σ) 卷积运算。

高斯函数：

尺度空间：

其中，*表示卷积计算。

3.2.3 高斯模糊

高斯模糊是在Adobe Photoshop等图像处理软件中广泛使用的处理效果，通常用它来减小图像噪声以及降低细节层次。这种模糊技术生成的图像的视觉效果是好像经过一个半透明的屏幕观察图像。

3.2.4 高斯金字塔

高斯金字塔的构建过程可分为两步：（1）对图像做高斯平滑；（2）对图像做降采样。

为了让尺度体现其连续性，在简单下采样的基础上加上了高斯滤波。一幅图像可以产生几组（octave）图像，一组图像包括几层（interval）图像。

高斯图像金字塔共o组、s层，则有：

σ——尺度空间坐标；s——sub-level层坐标；σ0——初始尺度；S——每组层数（一般为3~5）。

最后可将组内和组间尺度归为：

i——金字塔组数；n——每一组的层数。

3.3 检测尺度空间极值

3.3.1 DoG（Difference of Gaussian）函数

为了有效提取稳定的关键点，利用不同尺度的高斯差分核与卷积生成。

DOG在计算上只需相邻高斯平滑后图像相减，因此简化了计算！

3.3.2 DoG高斯差分金字塔

对应DOG算子，需构建DOG金字塔。可以通过高斯差分图像看出图像上的像素值变化情况。（如果没有变化，也就没有特征。特征必须是变化尽可能多的点。）DOG图像描绘的是目标的轮廓。

3.3.3 DOG局部极值检测

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

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

3.4 关键点的精确定位

以上方法检测到的极值点是离散空间的极值点，我们可以通过拟合三维二次函数来精确定位关键点的位置与尺度，由于DoG函数在图像边缘有较强的边缘响应，因此需要排除边缘响应，增强稳定性和抗噪声能力。

3.4.1 特征点的精确定位

子像素插值（Sub-pixel Interpolation）：利用已知的离散空间点插值得到的连续空间极值点的方法。

对尺度空间DOG函数进行曲线拟合，以提高关键点的稳定性，减小误差。

利用DoG函数在尺度空间的Taylor展开式(拟合函数)为：

其中 $X = (x, y, \sigma)^T$ ，求导并让方程等于零，可以得到极值点的偏移量为：

对应极值点，方程的值为：

其中 $X = (x, y, \sigma)^T$ 代表相对插值中心的偏移，当它在任意维度上的偏移量大于 0.5 时，意味着插值中心已经偏移到它的临近点上，所以必须改变当前关键点的位置。同时在新的位置上反复插值直到收敛；有也可能超出所设定的迭代次数或者超出图像边界的范围，此时这样的点应该删除。另外， $\left | D(x) \right |$ 过小的点易受噪声的干扰而变得不稳定，所以将 $\left | D(x) \right |$ 小于某个经验值的极值点删除。同时，在此过程中获取特征点的精确位置（原位置上加上拟合的偏移量）以及尺度。

3.4.2 剔除不稳定的边缘响应点

由于DoG函数在图像边缘有较强的边缘响应，因此需要排除边缘响应。 DoG函数的峰值点在边缘方向有较大的主曲率，而在垂直边缘的方向有较小的主曲率。主曲率可以通过计算在该点位置尺度的2×2的Hessian矩阵得到，导数由采样点相邻差来估计：

$H=\begin{bmatrix} D_xx& D_xy \\ D_xy & D_yy \end{bmatrix}$

其中 $D_{xx}, D_{xy}, D_{yy}$ ,是候选点邻域对应位置的差分求得的。

D的主曲率和H的特征值成正比。令 α ，β为特征值，则:

$\frac{Tr(H)^2}{Det(H)}=\frac{(\alpha +\beta )}{\alpha \beta }$

$Det(H)=D_x_x+D_y_y=\alpha+\beta$

$Tr(H)=D_x_xD_y_y-D_{xy}^2=\alpha \beta$

该值在两特征值相等时达最小。Lowe论文中建议阈值T为1.2，即:

$\frac{Tr(H)^2}{Det(H)}<T$

时保留关键点，反之剔除。

3.5 关键点方向分配

通过尺度不变性求极值点，可以使其具有缩放不变的性质。而利用键点邻域像素的梯度方向分布特性，可以为每个关键点指定方向参数方向，从而使描述子对图像旋转具有不变性。

通过求每个极值点的梯度来为极值点赋予方向。

像素点的梯度表示：

梯度幅值：

梯度方向：

计算得到梯度方向后，就是用直方图统计特征点邻域内像素对应的梯度方向和幅值。梯度方向的直方图的横轴是梯度方向的角度（梯度方向的范围是 0 到 360 度，直方图每 36 度一个柱，共 10 个柱），纵轴是梯度方向对应梯度幅值的累加，在直方图的峰值就是特征点的主方向。

关键点的主方向与辅方向：

关键点主方向：极值点周围区域梯度直方图的主峰值也是特征点方向
关键点辅方向：在梯度方向直方图中，当存在另一个相当于主峰值 80%能量的峰值时，则将这个方向认为是该关键点的辅方向。

这可以增强匹配的鲁棒性，Lowe的论文指出大概有15%关键点具有多方向，但这些点对匹配的稳定性至为关键。

3.6 关键点描述

下图是一个SIFT描述子事例。其中描述子由2×2×8维向量表征，也即是 2×2个8方向的方向直方图组成。左图的种子点由8×8单元组成。每一个小格都代表了特征点邻域所在的尺度空间的一个像素，箭头方向代表了像素梯度方向，箭头长度代表该像素的幅值。然后在4×4的窗口内计算8个方向的梯度方向直方图。绘制每个梯度方向的累加可形成一个种子点，如右图所示：一个特征点由4个种子点的信息所组成。

与求主方向不同，此时每个种子区域的梯度直方图在0-360之间划分为8个方向区间，每个区间为45度，即每个种子点有8个方向的梯度强度信息。在实际的计算过程中，为了增强匹配的稳健性，Lowe建议对每个关键点使用 4×4共16个种子点来描述，这样一个关键点就可以产生128维的SIFT特征向量。

通过对特征点周围的像素进行分块，计算块内梯度直方图，生成具有独特性的向量，这个向量是该区域图像信息的一种抽象，具有唯一性。

3.7 关键点匹配

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

模板图中关键点描述子：

实时图中关键点描述子：

任意两描述子相似性度量：

要得到配对的关键点描述子， $d(R_i,S_i)$ 需满足：

关键点的匹配可以采用穷举法来完成，但是这样耗费的时间太多，一般都采用kd树(一个平衡二叉树)的数据结构来完成搜索。搜索的内容是以目标图像的关键点为基准，搜索与目标图像的特征点最邻近的原图像特征点和次邻近的原图像特征点。

3.8 编程实现

3.8.1 检测感兴趣点

# 检测感兴趣点
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("D:\python\win64vlfeat\sift.exe " + 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 = 'picture1.jpg'
im = array(Image.open(imname).convert('L'))
process_image(imname, 'picture1.sift')
l1, d1 = read_features_from_file('picture1.sift')

figure()
gray()
subplot(131)
plot_features(im, l1, circle=False)
title('SIFT Feature')
subplot(132)
plot_features(im, l1, circle=True)
title('SIFT Feature Scale')

# 检测harris角点
harrisim = compute_harris_response(im)

subplot(133)
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 Corner Point')
show()

输出结果：

3.8.2 两张图像的SIFT特征匹配

# SIFT特征匹配
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("D:\python\win64vlfeat\sift.exe " + 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


# 双向对称版本的match()
def match_twosided(desc1, desc2):

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

运行测试：

im1f = 'picture1.jpg'
im2f = 'picture1.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()

输出结果：

3.8.3 地理标记图像匹配

from PIL import Image
from numpy import *
import os
import pydot
from PCV.localdescriptors import sift

# 获取图像列表
def get_imlist(path):
    return [os.path.join(path, f) for f in os.listdir(path) if f.endswith('.png')]

download_path = "D:/python/workSpace/ComputerVision/Test2/Pictures"
path = "D:/python/workSpace/ComputerVision/Test2/Pictures"

# 获取文件名列表
imlist = get_imlist(download_path)
nbr_images = len(imlist)
print(nbr_images)

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

# 复制值
for i in range(nbr_images):
    for j in range(i + 1, nbr_images):  # no need to
        matchscores[j, i] = matchscores[i, j]

# 可视化
# min number of matches needed to craete link
threshold = 2

# don't want the default directed graph
g = pydot.Dot(graph_type='graph')

for i in range(nbr_images):
    for j in range(i + 1, nbr_images):
        if matchscores[i, j] > threshold:
            # 图像对中的第一幅图像
            im = Image.open(imlist[i])
            im.thumbnail((100, 100))
            filename = path + str(i) + '.png'
            im.save(filename)  # 需要一定大小的临时文件
            g.add_node(pydot.Node(str(i), fontcolor='transparent',shape='rectangle', image=filename))

            # 图像对中的第二幅图像
            im = Image.open(imlist[j])
            im.thumbnail((100, 100))
            filename = path + str(j) + '.png'
            print(filename);
            im.save(filename)  # 需要一定大小的临时文件
            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('result.png')

输出结果：

3.8.4 实验过程中遇到的问题

问题一：

UserWarning: loadtxt: Empty input file: "out_sift.1.txt" f = loadtxt(filename)

IndexError: too many indices for array: array is 1-dinmensional，but 2 were indexed

—直报错类似的问题，还没解决!!!

错误原因：生成的sift文件为空，即特征提取为空。

解决办法：vlfeat0.9.20版本，下载了vcomp100.dll文件，下载链接https://cn.dll-files.com/vcomp100.dll.html，64位系统对应64位。

参考博客：vlfeat0.9.21提取sift特征为空_lrwin_bian的博客-CSDN博客

vlfeat安装参考教程：Python计算机视觉编程 - 安装

问题二：

关于Graphviz安装教程的一些问题：

第一步：确认在Graphviz官网下载并安装了对应操作系统（Windows，Linux，Mac）的Graphviz软件。官网地址：https://www.graphviz.org/。安装过程中有一个选项，建议选择添加至所有人，如下图：

第二步：将软件安装目录的bin文件添加到系统环境中。如我的电脑是Windows系统，Graphviz安装在D:\Program Files\Graphviz2.38，该目录下有bin文件，将该路径添加到系统环境变量中，即D:\Program Files\Graphviz2.38\bin。如果第一步选择了PATH添加至所有人，则这一步可以直接跳过！
第三步：使用pip命令安装以下模块：
                pip install pydot-ng
                pip install graphviz
                pip install pydot

注意安装顺序，一定要先安装grahiviz，再pip install pydot！！！
第四步：进入windows命令行界面，输入dot -version，然后按回车，如果显示graphviz的相关版本信息，则安装配置成功。如图：

依照上述步骤基本不会有问题，若仍报错"dot" not found in path，那么在提示目录下(如前文报错内容加粗目录)找到pydot.py文件的源码，在pydot.py中找到类Dot的定义，进行如下修改：

保存文件，如果运行仍报错"dot" not found in path，那就重启电脑！！！！

到这里，再运行就没出问题了。

DAMO开发者矩阵

DAMO开发者矩阵，由阿里巴巴达摩院和中国互联网协会联合发起，致力于探讨最前沿的技术趋势与应用成果，搭建高质量的交流与分享平台，推动技术创新与产业应用链接，围绕“人工智能与新型计算”构建开放共享的开发者生态。

更多推荐

DreamZero技术解析：当世界模型成为机器人“物理大脑”

DAMO开发者矩阵

硬实时·强生态：鸿道Intewell硬实时操作系统

DAMO开发者矩阵

【随手记】Covariant’s Brain Service和gRPC谷歌远程过程调用

高性能：HTTP/2 + Protobuf，高吞吐量、低延迟实时双向通信：机器人和 Brain Service 可持续交换数据支持流式数据：感知数据和动作指令都是连续流跨语言系统支持：简化分布式开发正是 gRPC 的这些特性，使得 Covariant 的机器人系统能够像“智能团队”一样高效协作，实现真正的工业自动化智能。