内容

1. 特征提取和匹配

1.1 背景知识

1.2 特征匹配基本流程

1.3 局部特征描述子

2. Harris角点检测

2.1 角点（corner points）

2.2 HARRIS角点检测基本思想

2.3 HARRIS检测：数学表达

2.4 角点响应函数

2.5 编程实现

2.5.1 角点检测代码实现

2.5.2 角点检测数据测试

3. SIFT特征匹配算法

3.1 SIFT算法实现步骤

3.2 关键点检测的相关概念

3.2.1 关键点（特征点）

3.2.2 尺度空间

3.2.3 高斯模糊

3.2.4 高斯金字塔

3.3 检测尺度空间极值

3.3.1 DoG（Difference of Gaussian）函数

3.3.2 DoG高斯差分金字塔

3.3.3 DOG局部极值检测

3.4 关键点的精确定位

3.4.1 特征点的精确定位

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

3.5 关键点方向分配

3.6 关键点描述

3.7 关键点匹配

3.8 编程实现

3.8.1 检测感兴趣点

3.8.2 两张图像的SIFT特征匹配

3.8.3 地理标记图像匹配

3.8.4 实验过程中遇到的问题

1. 特征提取和匹配

1.1 背景知识

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

关键点，也称为兴趣点，是在纹理中表示的点。关键点通常是对象边界方向突然改变的点或两个或多个边缘段之间的交点。它具有明确定义的位置或位于图像空间中的良好位置。即使图像域存在光照和亮度变化等局部或全局扰动，关键点仍然稳定，可以重复可靠地计算。除此之外，它还可以提供有效的检测。

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)$ ，我们可以使用二元泰勒近似来简化表达式，泰勒在 $(x,y)$ 处展开 $I(x+u,y+v)$ 函数得到：

$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，那就重启电脑！！！！

至此，再次运行就没有问题了。

文章出处登录后可见！

已经登录？立即刷新

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

1. 特征提取和匹配

1.1 背景知识

1.2 特征匹配基本流程

1.3 局部特征描述子

2.Harris角点检测

2.1 角点（corner points）

2.2 HARRIS角点检测基本思想

2.3 HARRIS检测：数学表达

2.4 角点响应函数

2.5 编程实现

2.5.1 角点检测代码实现

2.5.2 角点检测数据测试

3. SIFT特征匹配算法

3.1 SIFT算法实现步骤

3.2 关键点检测的相关概念

3.2.1 关键点（特征点）

3.2.2 尺度空间

3.2.3 高斯模糊

3.2.4 高斯金字塔

3.3 检测尺度空间极值

3.3.1 DoG（Difference of Gaussian）函数

3.3.2 DoG高斯差分金字塔

3.3.3 DOG局部极值检测

3.4 关键点的精确定位

3.4.1 特征点的精确定位

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

3.5 关键点方向分配

3.6 关键点描述

3.7 关键点匹配

3.8 编程实现

3.8.1 检测感兴趣点

3.8.2 两张图像的SIFT特征匹配

3.8.3 地理标记图像匹配

3.8.4 实验过程中遇到的问题

相关推荐