OpenCV入门学习笔记之Harris角点检测与SIFT特征匹配算法

xiaoxingxing

3 年前

1. 写在前面

这篇文章整理两个图像处理中非常重要的算法，一个是Harris角点检测算法，另一个是SIFT特征匹配算法，这两个算法本质上还是去找图像里面的关键特征点，帮助我们后续更好的理解图像以及做各种各样的分析。由于这两个算法涉及到的数学原理会比较多，而我刚入门，所以只是从使用的角度，简单的描述到底在做什么事情，至于详细的数学细节或者推导，这里不过多整理，以掉包能完成任务为首要目的啦。

首先，先介绍Harris角点检测算法，角点在图像中是很重要的特征，信息含量很高，那么如何找到一个图像里面的角点呢？这个算法就能轻松解决。但是呢？这个算法只考虑了旋转不变性，即一个角点旋转之后还是角点，没有考虑到尺度不变，尺度变化可能会导致角点变成边缘，所以有没有一种方法同时考虑这两个特性去找特征点呢？这就是SIFT算法，这个算法也是用来侦测和描述影像中的局部性特征，基于位置，尺度，旋转不变性在空间尺度中寻找特征点。这两个算法的使用场景非常广泛，比如图像配准，目标识别跟踪，图像对齐，图像拼接(全景图)，相机标定，三维场景重建等。最后，通过一个全景图拼接的Demo感受下这两个算法的魅力 😉

大纲如下：

Harris角点检测算法
SIFT特征匹配算法
特征匹配
案例Demo: 全景图拼接

Ok, let’s go!

2. Harris角点检测算法

角点特征能在保留图像重要特征的同时，有效减少信息数据量，算是图像中较好的特征，比边缘特征更好的用于定位。在图像处理中，检测角点特征的算法很多，最常用，最基础的就是Harris角点检测算法。

在说这个算法之前，我们先来感受一下什么是角点：

如果我们用人眼看，很容易理解，也就是图像中物体的一角。然后计算机可以看到，要定位这些角点并不容易。想让电脑看到，就必须明白一件事：这样的角点、边界点，和普通的屏幕点有什么区别，也就是图中红框、蓝框、黑框内的点。上图在数值上会有什么不同？

Harris角点检测认为：特征点具有局部差异性。如何描述呢？每一个点为中心，取一个窗口，比如窗口大小或者，那么窗口描述了这个特征点周围的环境。

如果特定窗口在图像的各个方向移动时，窗口中图像的灰度没有变化，那么窗口中不存在角点，比如蓝框
如果窗口向某个方向移动，窗口中图像的灰度变化很大，而其他方向没有变化，那么窗口中的图像可能是一条直线边界，比如黑框，而垂直运动剧烈变化，水平方向不变
如果这个特征的小窗口向四面八方移动，窗口内区域的灰度变化很大，则认为窗口存在角点。

说是这么说，但具体应该怎么找到角点呢？ Harris角点检测算法主要是三步：

当窗口同时向和两个方向移动时，计算窗口内像素值的变化
对于每个窗口，计算对应的角点响应函数
对于这个函数的阈值处理，如果，表示窗口对应一个角点特征

下面详细展开。

2.1 窗口移动，计算像素值变化量

这里的核心问题：如何确定哪些窗口导致灰度值变化较大？

假设像素位置在当前窗口的中心，该点的像素值为，如果窗口分别向和方向移动和步，到达新的位置，则像素值为，则是中心点移动所引起的灰度值变化。

那我们这是个窗口呀，比如，那就有9个点，窗口一移动的话，这9个点都有引起的灰度值的变化，而把这个都加和，就得到了窗口在各个方向上移动造成的像素灰度值变化，这个应该很好理解，公式如下：

这里有一个函数，表示窗口中每个像素的权重。可以设置为以窗口中心为原点的高斯分布。如果窗口中心的像素是角点，那么在窗口移动前后，中心点的灰度值变化非常强烈，则权重较大，说明该点对灰度变化的贡献很大。对于离窗口中心越远的点，灰度变化越小，因此权重越小，对灰度变化的贡献就越小。

这个公式是我们的目标函数。如果是角点，这个函数的值会比较大，所以我们最大化这个函数，得到图像中的角点。

But，上面这个函数计算会非常慢，比较涉及到了窗口内所有像素点的计算，所以，我们用泰勒，先对像素值函数进行近似。

其中，和是的偏微分，是图像中和方向的梯度图，可以通过cv2.Sobel()得到

将上式代入化简：

这由一个矩阵表示，即线生成的二次变换：

其中，矩阵如下：

所以最终的目标函数变为：

在这一点上，应该很好理解。无非就是直线生成中的泰勒近似和二次化简。注意在这里是一个协方差矩阵。主对角线可以看成是梯度在自身方向的方差，副对角线是与其他方向梯度的协方差，这是一个对称矩阵。

以下几点不是很好理解：二次函数本质上是一个椭圆函数，椭圆方程为：

可视化如下：

这里是实对称矩阵的特征值。这里可能不太好理解。我将尝试从基础转换的角度来解释它。首先，是一个实对称矩阵，那么它必须满足：

是的特征向量组合，是的特征值。将此公式代入上式：

这就变成了一个标准的椭圆方程，的特征值的平方根就代表了椭圆的长轴和短轴。那么特征值是什么？其实就是将原来的像素点映射到新的空间中（映射规则为），每组基方向的梯度方差或变化程度。所以我们可以用这个特征值来衡量像素在某个方向上的波动程度。这里如果不明白，可以看我的文章，补一下向量表示和基变换的知识。

因此，通过以上操作，将某个方向的灰度变化程度转化为矩阵的特征值。

2.2 窗口的角点响应函数

这里定义了每个窗口的角响应函数

这里是一个介于之间的经验常数。

2.3 非极大值抑制

根据 R 的值，将这个窗口所在的区域划分为平面、边缘或角点。为了得到最优的角点，我们还可以使用非极大值抑制。

当然，既然特征值和决定了的值，我们可以用特征值来判断一个窗口是平面、边还是角：

平面：窗口在平坦的区域滑动，窗口内的灰度值基本没有变化，所以很小，水平和垂直方向的变化比较小，即和小，和小。
边：，仅在水平或垂直方向有较大变化，和有较大，即或
角点：值大，水平和垂直方向变化大的点，即和大，和大。

一张图片胜过千言万语：

Harris 角点检测的结果是带有这些分数的灰度图像，设定一个阈值，分数大于这个阈值的像素就对应角点。

注意：Harris 检测器具有旋转不变性，但不具有尺度不变性，也就是说尺度变化可能会导致角点变为边缘，如下图所示：

So，如何找到旋转以及尺度不变的特征点呢？这就是SIFT算法干的事情了，但是介绍之前，先看看OpenCV中角点检测算法咋用。理论一大推，但是用起来一个函数搞定。

2.4 OpenCV中的角点检测

这里的函数是cv2.cornerHarris()：

img: 数据类型为float32的输入图像
blockSize: 角点检测中指定区域的大小，即窗口大小
ksize: Sobel求导中使用的窗口大小，需要用sobel算子求梯度，这里是设置sobel算子求导的窗口
k：取值参数为[0.04, 0.06]

看一个例子：

img = cv2.imread('img/test_img.jpg')
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# 角点检测
dst = cv2.cornerHarris(gray, 2, 3, 0.04)  # 这个是每个像素点的E值，即平移后灰度级变换程度值

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

以下是角点检测结果：

3. SIFT算法

Scale Invariant Feature Transform(SIFT): 尺度不变特征转换用来侦测与描述影像中的局部性特征，基于位置，尺度和旋转不变性在空间尺度中寻找极值点。

特征：

SIFT特征是图像的局部特征，其对旋转、尺度缩放、亮度变化保持不变性，对视角变化、仿射变换、噪声也保持一定程度的稳定性
唯一性好，信息丰富，适用于数据库中海量特征的快速准确匹配
多量性，即使少数几个物体也可以产生大量的SIFT特征向量
高速和可扩展性，特征可以很容易地与其他向量结合

解决问题：目标自身状态，场景所处环境和成像器材的成像特性等因素影响图像配准/目标识别跟踪的性能。SIFT算法一定程度可解决：

目标旋转、缩放、平移
图像仿射/投影变换
灯光效果
目标遮挡
杂物场景
噪音

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

SIFT算法主要下面四步：

尺度空间极值检测：在所有尺度上搜索图像位置，并通过高斯微分函数识别潜在的尺度和旋转不变特征点
关键点定位：在每个候选位置，位置和规模由一个拟合良好的模型确定。关键点是根据它们的稳定性来选择的。
方向确定：根据图像的局部梯度方向为每个关键点位置分配一个或多个方向。对图像数据的所有后续操作都相对于关键点的方向、比例和位置进行转换，从而为这些转换提供不变性。
关键点描述：在每个关键点周围的邻域中，测量所选尺度下图像局部性的梯度。这些梯度被转换成允许相对较大的局部形状变形和光照变化的表示。

由于这个算法理论上稍微复杂一些，下面是一个简单的安排。首先，我们要先了解这个算法是做什么的。

我明白了：找到图像中那些具有旋转、平移和尺度不变性的特征点，最后将它们表示为一个向量。

它究竟是如何完成的？大致内容如下。

3.1 图像的尺度空间

在一定范围内，无论物体大小，人眼都能分辨。但是，计算机很难拥有同样的能力。为了让机器对不同尺度的物体有统一的认知，需要考虑不同尺度图像中存在的特征。

尺度空间的获取通常使用高斯模糊来实现。

高斯模糊是一种图像滤波器，它利用正态分布（高斯函数）计算出一个模糊模板，并利用该模板与原始图像进行卷积运算，达到对图像进行模糊处理的目的。公式如下：

参数为标准差。指定值越大，像素的变化越大，与原始图像的偏差越大，即图像越模糊。模糊半径，指模板元素到模板中心的距离，如果二维模板的尺寸维度为，则模板上元素对应的高斯计算公式：

由具有非零分布的像素组成的卷积矩阵与原始图像进行变换。每个像素的值是周围相邻像素值的加权平均值。原始像素的值具有最大的高斯分布值，因此具有最大的权重，并且随着与原始像素的距离越来越远，相邻像素的权重越来越小。这种模糊比其他均衡模糊过滤器更好地保留了边缘效果。

不同的模糊效果如下：

3.2 多分辨率金字塔

尺度空间使用高斯金字塔表示，尺度规范化LoG(Laplacion of Gaussian)算子具有真正尺度不变性，Lowe使用高斯差分金字塔近似LoG算子，在尺度空间检测稳定关键点。

实现尺度空间时，用高斯金字塔表示。高斯金字塔的构建主要分为两部分：

在图像上做不同尺度的高斯模糊
对图像进行下采样以构建金字塔

高斯金字塔应该很熟悉：

图像的金字塔模型是指对原始图像进行连续下采样得到一系列不同大小的图像，从大到小，从下到上，形成的塔形模型。原始图像是金字塔的第一层，每次下采样得到的新图像是金字塔的一层（每层一张图像），每个金字塔共有层。金字塔的层数由图像的原始大小和塔顶图像的大小决定。计算公式如下：

其中，是原始图像的大小，是塔顶图像最小维度的对数。

3.3 DoG空间极值点检测

在实际计算中，将高斯金字塔中每组图像的相邻上下层相减，得到高斯差分图像。

这个就是高斯金字塔中的每一层的图片，进行差分，这样得到的结果中，像素点相差较大的位置，就是不同尺度的图片之间的不同。DoG公式定义如下：

这个公式也很容易理解。差分结果等于第一个高斯滤波器的结果与第二个高斯滤波器的结果之差。

为了寻找尺度空间的极值点，每个像素点要和其图像域(同一尺度空间)和尺度域(相邻的尺度空间)的所有相邻点进行比较，当其大于(或者小于)所有相邻点时，该点就是极值点。如下图所示，中间的检测点要和其所在图像的邻域的8个像素点，以及其相邻的上下两层的领域的18个像素点，共26个像素点进行比较。

由于要在相邻尺度进行比较，如上上面那个图，每组含4层高斯差分金字塔，只能中间两层中进行两个尺度的极值点检测，为了在每组中检测个尺度的极值点，则DOG金字塔每组需要层图像，而DOG金字塔由高斯金字塔相邻两层相减得到，则高斯金字塔每组需层图像，实际计算时在3-5之间。

当然这样产生的极值点并不全都是稳定的特征点，因为某些极值点响应较弱，而且DoG算子会产生较强的边缘响应。

值得一提的是，选取的高斯差分金字塔极值点只是候选特征点。虽然高斯差分金字塔的极值点已经能够更好地表示图像的特征并且具有尺度不变性，但是在选择过程中没有考虑图像特征点对图像噪声的鲁棒性，确定的图像特征点在这种方式在实际应用中使用。容易出现图像匹配不当等问题。

3.4 关键点的精确定位

以上方法检测到的极值点是离散空间的极值点，以下通过拟合三维二次函数来精确确定关键点的位置和尺度，同时去除低对比度的关键点和不稳定的边缘响应点(因为DoG算子会产生较强的边缘响应)，以增强匹配稳定性、提高抗噪声能力。

为了更好的理解一点，我们先来看看如何在一维的离散值拟合一条曲线：

将其置于 3D 中：

这里的表示相对窗口中心的偏移量，等价于上一节中的。如果对上面进行求导，并令导数等于0，就可以得到较为准确的极值点和极值。

这里将每个候选极值点求出相应的导数，然后代入，得到，把结果值非常小的(比如小于0.03)的先进行首轮剔除。

3.5. 消除边界响应

DoG算子会产生较强的边缘响应，即容易保留边界点，所以下面需要剔除不稳定的边缘响应点。具体做法如下：

首先，获取每个特征点出的Hessian矩阵，这里其实和角点检测那里是一样的：

设为最大特征值，为最小特征值，求下界值：

通过相邻采样点之间的差异来估计导数。比值越大，两个特征值的比值就越大，即一个方向的梯度值越大，而另一个方向的梯度值越小，边缘正是如此。因此，为了消除边缘响应点，该比率需要小于某个阈值。因此，对于每个特征点，传递以下公式：

建立这种不等式的关键点被保留，反之亦然。论文

3.6 特征点的主方向

根据上述操作，可以找到一些较好的候选关键点，但是在构建高斯差分金字塔选择图像特征点时，算法考虑了关键点的尺度不变性。对于图像特征，旋转不变性与尺度不变性同等重要。

为了使描述旋转不变，需要利用图像的局部特征为每个关键点分配一个主方向。保持与旋转前相同的特征。

该算法计算关键点指定大小的场中所有点的梯度方向和赋值，并将所有梯度方向对应的赋值和计算为关键点周围邻域梯度方向的直方图。梯度的模值和方向计算如下：

这样，对于每个关键点，都有三条信息，即位置、尺度和方向。但关键点的主要方向是什么？

对于一个关键点，计算其邻域内每个点的梯度直方图，

下图直方图为简化版本（只有8个bin，实际操作时算法会统计从0到360°步长为10°的共计36个梯度方向的幅值和，共有36个bin）。梯度方向直方图中最高的bin对应的方向即定义为该关键点的主方向，若存在任一方向的幅值大于主方向幅值的80%，则将其作为辅方向。所以，对于同一个关键点，可能有多个方向，这种情况在相同位置和尺度将会有多个关键点被创建但方向不同。实际编程实现中，就是把该关键点复制成多份关键点，并将方向值分别赋给这些复制后的关键点。

直方图的横轴是方向，纵轴是邻域内每个点对应的方向梯度的累积和。

3.7 生成特征描述

在得到特征点的二维位置、尺度位置、主方向的具体信息后，算法需要解决的最后一个问题就是生成关键点信息的描述符，即用一个向量来描述图像中的特征点信息。它不随各种变化而变化，如光照变化、视角变化等。该描述符不仅包括关键点，还包括对其有贡献的关键点周围的像素，并且描述符应具有较高的唯一性以提高特征点正确匹配的概率。

为了保证特征向量的旋转不变性，需要以特征点为中心，在附近场内将坐标轴旋转角，即以坐标轴为主方向旋转的特征点。

每个像素的坐标变换使用以下公式：

旋转完之后，算法将特征点周围的邻域分为4个的区域，再将的区域划为区域，即每个小区域为的范围。统计每个区域的梯度方向直方图（直方图为8个bin，8个方向），故共计个bin。对应生成128维向量（值为梯度方向的幅值）。该128维向量即该点的特征描述子。

每个像素的梯度幅值和方向(箭头方向代表梯度方向，长度代表梯度幅值)，然后利用高斯窗口对其加权运算(离中心点距离不一样)，最后在每个的小块上绘制8个方向的梯度直方图，计算每个梯度方向的累加值，即可形成一个种子点，即每个特征点的由个种子点组成，每个种子点由8个方向的向量信息。

论文中建议对每个关键点使用共16个种子点来描述，每个种子点由个方向的梯度直方图描述，即8维向量，这样一个关键点就会产生128维的SIFT特征向量

OK，到这里，就把SIFT算法的大致流程整理了下，不过上面有些乱，下面集中总结下：

对于原始图像，先用高斯滤波做模糊操作的尺度变换，通过改变高斯核的标准差，得到图像
为了让后序特征表达更丰富，把这张图做成了高斯金字塔
接下来，对于高斯金字塔的每一个level的张图片，相邻图片进行差分运算，得到高斯差分图像，张
基于高斯差分图像，找到候选极值点，所谓极值点就是比较差分图像中的每个像素点，而不是图像的第一层和最后一层，比较不同的尺度空间（上层和上层）下相邻两层邻域）和周围像素在同一尺度空间的值和该像素的当前值，得到最大值点作为候选关键点。
定位精确的关键点，对尺度空间中DOG函数进行曲线拟合，计算其极值点，从而实现关键点的精确定位。这里用到了泰勒近似，拟合出函数来之后，把当前关键点代入，得到的函数值过小的(比如小于0.03)的关键点去掉。
消除边缘效应，得到剩余关键点在特征点处的梯度矩阵，用有限差分法求导。得到这个矩阵后，根据特征值的比值进行过滤。水平梯度类似于垂直梯度。如果保留，则差异很大，说明是一条边。这个关键点应该去掉，也就是消除边缘效应。
通过以上步骤，选择了关键点，但只是从尺度不变性的角度出发，再考虑平移不变性，为每个关键点选择一个主方向
这个就是考虑关键点邻域内的所有像素点，用一个直方图去统计这个邻域内所有像素点的梯度方向以及梯度累加值，梯度方向由于是360度，这里进行了分桶操作，划分成了8个bins。拿到直方图之后，把梯度累加值最大的那个方向作为关键点的主方向，当然次大的还可以作为辅方向，这样每个关键点就同时有了位置，尺度，方向三个特征描述(如果关键点还有辅方向的，就需要把这个观测点复制一份，保持位置，尺度不变，该变方向即可)。
每个观察点都有三个属性描述。接下来需要转化为特征，也就是用一个向量来描述
首先，对于每个观测点值，首先修改邻域内所有像素点的坐标值。修改的原则是坐标轴要与当前观测点的主方向一致。
修改完毕之后，统计这个邻域内的像素点的梯度直方图，依然是8个bins，每个bins中是像素点梯度的累加值，这样对于一个关键点，锁定一个邻域就能得到8个值。这算是一个种子区间内的特征描述
我们实际做的时候，首先给一个大的邻域一个关键点，然后把大的邻域划分成小的邻域作为种子区间，然后得到种子区间的特征。
关键点最后的特征描述，就是所有种子矩阵特征的Concat值，上面是一个128维的向量表示，这样就得到了关键点尺度和平移不变的特征描述
金字塔中的所有关键点获取特征描述并返回结果

你可以看到下面的代码。

3.8 OpenCV下的SIFT算法

对于SIFT算法， OpenCV中直接也是一个函数搞定。

img = cv2.imread('img/test_1.jpg')
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# opencv版本高于3.4.3， 这个sift算法使用变成cv2.SIFT_create()， 在这之前的版本是cv2.xfeatures2d.SIFT_create()

# SIFT检测器
sift = cv2.SIFT_create()

# 找出图像中的关键点
kp = sift.detect(gray, None)

# 在图中画出关键点
img = cv2.drawKeypoints(gray, kp, img)  

# 计算关键点对应的SIFT特征向量
kp, des = sift.compute(gray, kp)   

# 上面的两步，也可以用下面的一步
sift = cv2.SIFT_create()
kp1, des1 = sift.detectAndCompute(gray, None)   # 这里还能一步到位，直接算出关键点以及关键向量

看看效果：

4. 特征匹配

得到图像中的关键点后，关键点的特征也可以用向量来描述，给出两张图像后，如何分辨出两张图像中哪些关键点相似？这是看关键点向量之间的差异，也是特征匹配在做的事情。

4.1 Brute-Force蛮力匹配

暴力匹配，两个图像中的关键点的特征向量，一个个的计算差异，即两层for循环了。这里直接看怎么用：

首先，读入两张图像，然后得到它们各自的关键点和关键点向量。

img1 = cv2.imread('img/box.png', 0)
img2 = cv2.imread('img/box_in_scene.png', 0)

# sift算法，得到每张图的关键点以及关键向量
sift = cv2.SIFT_create()
kp1, des1 = sift.detectAndCompute(img1, None)   # 这里还能一步到位，直接算出关键点以及关键向量
kp2, des2 = sift.detectAndCompute(img2, None)

特征匹配执行如下：

1对1的匹配：即图片A中的一个向量匹配图片B中的一个向量# crossCheck表示两个特征点要互相匹配，例如A中的第i个特征点与B中的第j个特征点最近的，并且B中的第j个特征点到A中的第i个特征点也是最近的# NORM_L2: 归一化数组的欧几里得距离，如果其他特征计算方法需要考虑不同的匹配计算方式bf = cv2.BFMatcher(crossCheck=True)# 1对1的匹配matches = bf.match(des1, des2)matches = sorted(matches, key=lambda x: x.distance)
如果这里指定属性crossCheck为False，得到的结果是604，如果为True，得到的结果是206，所以我基于这个，盲猜下暴力匹配以及这个属性的意义。
暴力匹配的话，就是对于A图像中的每个观测点的向量，我遍历一遍B图像中每个观测点的向量，然后求欧几里得距离，拿到最小的。这样对于A图像中每个观测点，就得到了B图像里面的最佳匹配。
然后再对于B图像中的每个观测点，也同样用上面的方式走一遍，这样对于B图像中的每个观测点，也得到了A图像中的最佳匹配。
此时，如果是：
bf.match(des1, des2): 返回的就是A图像中的每个观测点的最佳匹配，个数是A中关键点的个数
bf.match(des2, des1): 返回的是B图像中每个观测点的最佳匹配，个数是B中关键点的个数
如果设置crossCheck为True的话：
bf.match(des1, des2)和 bf.match(des2, des1) 都是260，表示的其实是A和B中，关键点互相匹配的那些，比如A中的第j个观测点，最近点是B中的i，那么B中的i，也要对应A中的j，即 VS 。 img3 = cv2.drawMatches(img1, kp1, img2, kp2, matches[:10], None, flags=2)
画10个最近的匹配看看：
K对最佳匹配: 对于一个关键点，找K个最相似的向量。bf = cv2.BFMatcher()matches = bf.knnMatch(des1, des2, k=2) # 相当于对于A中的每个关键点，不是找某一个最相似，而是K相似
这里可以根据近邻之间的相似比例对关键点筛选 good = []for m, n in matches: if m.distance < 0.75 * n.distance: good.append([m])
也可以可视化一下： img3 = cv2.drawMatchesKnn(img1, kp1, img2, kp2, good, None, flags=2)
结果如下：

这里会发现一个问题，对于上面的匹配结果，大部分结果还可以，但有某些匹配错误的点，这种情况怎么弥补呢？可以使用RANSAC(Random sample consensus)随机抽样一致算法。这个东西可以简单看下原理。

4.2 RANSAC(Random sample consensus)算法

该算法也是一种拟合算法，与最小二乘法对比如下：

思路如下：首先选择初始样本点进行拟合，给出一个容差范围，继续迭代。

每次拟合后，在公差范围内都有对应的数据点，找到数据点数最多的案例就是最终的拟合结果。因此，该算法在拟合时，先随机采样初始点，然后在拟合时考虑容差范围（可以拟合的数据点个数）。

那么，这东西干啥用呢？下面全景拼接的Demo会用到。

5. 全景拼接Demo

全景拼接想必大家都玩过。此功能在相机中可用。其原理大概如下。假设有两张图片要拼接，大概的流程如下：

通过SIFT算法拿到两张图片的关键点以及关键向量
基于关键向量的特征匹配
根据特征匹配点，对图片进行一些仿射变换，比如平移、选择、缩放等，这里是为了保证可以无缝连接。如果不加工，则不得缝合。
既然是仿射变换，本质就是将变换矩阵相乘，所以核心就是求解这个变换矩阵
如果需要做变换，目的就是求变换矩阵，至少有4对匹配好的特征点才行。
那么思路就是，先用随机采样算法采出四对匹配点，然后解8个方程得到上面的变换矩阵，然后基于这个变换矩阵算其他的匹配点是否能匹配上，定义一个loss函数为匹配上的点的对数，这样就能基于当前的四对匹配点算出一个loss值。然后再随机采样四对，再求变换矩阵，再看看其他匹配点能匹配上的个数，得到损失。依次迭代多次，取loss函数最小的那个变换矩阵。
找到变换矩阵，先变换一张图片，再拼接。

这么说比较抽象，下面通过一个Demo把这个过程带起来。

5.1 读入图像，拿到SIFT特征

代码显示如下：

image_left = cv2.imread('img/left_01.png')
image_right = cv2.imread('img/right_01.png')

两张图片如下：

显然，这种拼接是不能拼接的。

获取下面的关键点和特征向量

def detectAndDescribe(image):
    # 将彩色图片转成灰度图
    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
    
    # SIFT生成器
    destriptor = cv2.SIFT_create()
    (kps, features) = destriptor.detectAndCompute(image, None)
    
    # 结果转成numpy数组
    kps = np.float32([kp.pt for kp in kps])
    return (kps, features)
 
# 检测A， B图片的SIFT特征关键点，得到关键点的表示向量
(kps_left, features_left) = detectAndDescribe(image_left)   # kpsA (关键点个数， 坐标)  features(关键点个数，向量)
(kps_right, features_right) = detectAndDescribe(image_right)

5.2 特征匹配

这里还是要写一个函数：

def matchKeyPoints(kpsA, kpsB, featuresA, featuresB, ratio=0.75, reprojThresh=4.0):
    # 建立暴力匹配器
    matcher = cv2.BFMatcher()
    
    # KNN检测来自两张图片的SIFT特征匹配对
    rawMatches = matcher.knnMatch(featuresA, featuresB, 2)
    
    matches = []
    for m in rawMatches:
        # 当最近距离跟次近距离的比值小于ratio时，保留此配对
        # (<DMatch 000001B1D6B605F0>, <DMatch 000001B1D6B60950>) 表示对于featuresA中每个观测点，得到的最近的来自B中的两个关键点向量
        if len(m) == 2 and m[0].distance < m[1].distance * ratio:   
            # 存储两个点在featuresA, featuresB中的索引值
            matches.append([m[0].trainIdx, m[0].queryIdx])    # 这里怎么感觉只用了m[0]也就是最近的那个向量啊，应该没用到次向量
            # 这个m[0].trainIdx表示的时该向量在B中的索引位置， m[0].queryIdx表示的时A中的当前关键点的向量索引
    
    # 当筛选后的匹配对大于4时，可以拿来计算视角变换矩阵
    if len(matches) > 4:
        # 获取匹配对的点坐标
        ptsA = np.float32([kpsA[i] for (_, i) in matches])
        ptsB = np.float32([kpsB[i] for (i, _) in matches])
        
        # 计算视角变换矩阵  这里就是采样，然后解方程得到变换矩阵的过程
        (H, status) = cv2.findHomography(ptsA, ptsB, cv2.RANSAC, reprojThresh)
        
        return (matches, H, status)
    
    # 匹配结果小于4时，返回None
    return None

主要逻辑是从图片B中给图片A中的关键点拿最近的K个匹配向量，然后基于规则筛选，保存好匹配好的关键点的两个索引值，通过索引值取到匹配点的坐标值，有了多于4对的坐标值，就能得到透视变换矩阵。这里返回的主要就是那个变换矩阵。

# 匹配两张图片的所有特征点，返回匹配结果 注意，这里是变换right这张图像，所以应该是从left找与right中匹配的点，然后去计算right的变换矩阵
M = matchKeyPoints(kps_right, kps_left, features_right, features_left)
if not M:
    # 提取匹配结果
    (matches, H, status) = M

这里要一定要注意好，到底是对哪张图片做变换，比如给图片B做变换，那么就从A中找与B中特征点匹配的特征向量，求的是让图片B变换的透视矩阵。这里是对right做变换，所以从left中给right的关键点找匹配点，给right计算透视矩阵，接下来，变换right

# 图片right进行视角变换， result是变换后的图片
result = cv2.warpPerspective(image_right, H, (image_left.shape[1] + image_right.shape[1], image_right.shape[0]))

cv_imshow('result', result)

下面看一下改造后的效果：

接下来就可以进行拼接了：

# 将图片A传入result图片最左端
result[0:image_right.shape[0], 0:image_right.shape[1]] = image_left

结果如下：

这样就拼接在一起了，除了亮度不一样，效果还是不错的。

5.3 画出匹配图像

这里也可以在图像上绘制匹配向量：

def drawMatches(imageA, imageB, kpsA, kpsB, matches, status):
    # 初始化可视化图片，将A、B图左右连接到一起
    (hA, wA) = imageA.shape[:2]
    (hB, wB) = imageB.shape[:2]
    vis = np.zeros((max(hA, hB), wA + wB, 3), dtype="uint8")
    vis[0:hA, 0:wA] = imageA
    vis[0:hB, wA:] = imageB

    # 联合遍历，画出匹配对
    for ((trainIdx, queryIdx), s) in zip(matches, status):
        # 当点对匹配成功时，画到可视化图上
        if s == 1:
            # 画出匹配对
            ptA = (int(kpsA[queryIdx][0]), int(kpsA[queryIdx][1]))
            ptB = (int(kpsB[trainIdx][0]) + wA, int(kpsB[trainIdx][1]))
            cv2.line(vis, ptA, ptB, (0, 255, 0), 1)

    # 返回可视化结果
    return vis

vis = drawMatches(image_left, image_right, kps_left, kps_right, matches, status)

结果如下：

6. 小总

这篇文章主要是整理了在图像处理中重要且常用的找特征点的两个算法Harris和SIFT算法，包括算法的数学原理以及如何使用，然后是整理了下特征匹配与一致性采样算法，这俩东西其实为了后面的透视变换矩阵服务。最后通过一个全景图拼接的Demo感受了下算法的魅力。

当然，这只是冰山下的一个小角落，因为后面我们会在一个项目中用到这个知识来识别停车场的车位。另外，这些方法还是通用的算法，所以我也想静下心来待会儿复习。

参考：

OpenCV入门到实战课程
角点检测：Harris 与 Shi-Tomasi
SIFT算法详解
SIFT算法简述及Python标记SIFT特征检测实践

文章出处登录后可见！

已经登录？立即刷新