SIFT 图像特征匹配算法

SIFT简介

SIFT（Scale-Invariant Feature Transform，尺度不变特征变换匹配算法）是一种高效区域检测算法。

SIFT算法可以解决的问题：

1. RST（Rotation Scale Translation）：图像的旋转、缩放、平移等变换
2. 图像的仿射变换（Affine Transformation）：图像的拉伸、压缩、扭曲等变换
3. 图像的光照变化
4. 目标部分遮挡
5. 杂物和噪声干扰

SIFT算法原理

1. 高斯滤波

一个图像的尺度空间$L$是高斯函数$G(x,y,\sigma )$与图像$I(x,y)$的卷积：

$$\begin{array}{rl}L(x,y,\sigma )& =G(x,y,\sigma )\ast I(x,y)\\ G(x,y,\sigma )& ={\displaystyle \frac{1}{2\pi {\sigma }^2}}{e}^{-\frac{(x-x_0{)}^2+(y-y_0{)}^2}{2{\sigma }^2}}\end{array}$$

其中，$G(x,y,\sigma )$是高斯函数，$\sigma$是尺度因子，$\ast$是卷积运算，$x_0,y_0$是高斯函数的中心（卷积核的中心）。

$\sigma$越小，图像被平滑的越少，图像的细节越多。
小尺度对应于图像的细节，大尺度对应于图像的整体特征。

2. 高斯金字塔

高斯金字塔包含多组（Octave）图像，每组图像由多层（Interval）不同尺度的高斯模糊图像组成。

构建过程：

1. 将原图像扩大到两倍，作为第0组的第0层
2. 对第0组的第0层进行${\sigma }_0$的高斯模糊，得到第0组的第1层
3. 选定一个比例系数$k$，对第0组的第1层进行$k{\sigma }_0$的高斯模糊，得到第0组的第2层
4. 对第0组的第2层进行$k^2{\sigma }_0$的高斯模糊，得到第0组的第3层，以此类推…
5. 将上一组倒数第3层图像做比例因子为2的降采样，得到下一组的第0层，重复步骤2-4

组数计算公式：$O=\lfloor {\log }_2(min(M,N))\rfloor -3$

• $M,N$：原图像的长宽

层数公式：$S=n+3$

• $S$是每组的层数（自设）
• $n$是最终想在差分金字塔中提取极值点的层数（见差分金字塔）

平滑因子公式：$\sigma (o,r)={\sigma }_0{2}^{o+\frac{r}{n}}$，其中$o$是组数，$r$是层数

• SIFT算法中，${\sigma }_0=1.6$。但由于相机具有初始模糊${\sigma }_0=0.5$，实际${\sigma }_0=\sqrt{{1.6}^2-{0.5}^2}=1.52$
• 第o层的初始平滑因子$\sigma (o,0)=2^o{\sigma }_0$
• $k=2^{1/n}$

金字塔内各图片的关系：

• 在同一组内，不同层图像的尺寸是一样的，后一层图像的高斯平滑因子σ是前一层图像平滑因子的k倍
• 在不同组内，后一组第一个图像是前一组倒数第三个图像的二分之一采样，图像大小是前一组的一半

3. 高斯差分金字塔与极值点检测

3.1 高斯差分金字塔

对高斯金字塔的每一组图像，计算相邻两层图像的差分，得到高斯差分金字塔(DoG, Difference of Gaussian)。

$$D(x,y,\sigma )=(G(x,y,k\sigma )-G(x,y,\sigma ))\ast I(x,y)$$

3.2 估计极值点

将每个像素点和它在3×3×3的邻域内的26个像素点进行比较。
如果这个点是这27个点中的最大值或最小值，则认为这个点是一个极值点。

DoG金字塔是离散的(因为尺度空间和像素点都是离散的)，所以找到的极值点不太准确的，很大可能在真正极值点附近，因此需要对每个极值点进行插值。

3.3 泰勒展开求精确极值点

对每个极值点$X_0(x_0,y_0,{\sigma }_0{)}^T$，进行三元二次泰勒展开：

$$\begin{array}{rl}f(X)& =f(x_0)+{{\displaystyle \frac{\partial f}{\partial X}}}^T\hat{X}+{\displaystyle \frac{1}{2}}{\hat{X}}^T{\displaystyle \frac{{\partial }^{2}f}{\partial X^2}}\hat{X}\\ \hat{X}& =X-X_0\\ X& =(x,y,\sigma {)}^T\end{array}$$

将$f(X)$对$X$求导：

$${\displaystyle \frac{\partial f(X)}{\partial X}}={{\displaystyle \frac{\partial f}{\partial X}}}^T+{\displaystyle \frac{1}{2}}({\displaystyle \frac{{\partial }^{2}f}{\partial X^2}}+{{\displaystyle \frac{{\partial }^{2}f}{\partial X^2}}}^T)\hat{X}={{\displaystyle \frac{\partial f}{\partial X}}}^T+{\displaystyle \frac{{\partial }^{2}f}{\partial X^2}}\hat{X}$$

导数为0时，解得极值点：

$$\hat{X}=-{{\displaystyle \frac{{\partial }^{2}f}{\partial X^2}}}^{-1}{\displaystyle \frac{\partial f}{\partial X}}$$

代入$f(X)$得到极值：

$$f(X)=f(x_0)+{\displaystyle \frac{1}{2}}{{\displaystyle \frac{\partial f}{\partial X}}}^T\hat{X}$$

3.4 去除低对比度点

舍去$|f(X_0)|<{\displaystyle \frac{T}{n}},T=0.04$的极值点，因为这些点对应的差分值太小，对比度低，不具有代表性。

3.5 去除边缘响应

DoG在边缘会有比较大的值（边缘响应强），但边缘不一定能够提供稳定的特征，因此需要去除边缘响应强的点。

对于每个极值点，计算Hessian矩阵的迹和行列式：

$$\begin{array}{rl}H& =\left[\begin{array}{ccc}{D}_{xx}& D_{xy}{D}_{xy}& D_{yy}\end{array}\right]\\ Tr(H)& =D_{xx}+D_{yy}=\alpha +\beta \\ Det(H)& =D_{xx}{D}_{yy}-(D_{xy}{)}^2=\alpha \beta \end{array}$$

其中：

• $\alpha$为较大的特征值，$\beta$为较小的特征值
• $\alpha =r\beta$，$r$是一个大于1的实数

过滤：

1. $Det(H)<0$
2. ${\displaystyle \frac{(Tr(H){)}^2}{Det(H)}}\ge (r+1{)}^2/r$，论文推荐值$r=10$

4. 关键点方向分配

1. 0~360度划分为36 bins，每个bin为10度
2. 在高斯金字塔找到关键点对应位置，以它为圆心，半径为$1.5\sigma$画圆
3. 统计圆中所有像素的梯度方向、梯度幅值（模），直方图平滑处理，用$1.5\sigma$进行高斯加权
   • 平滑处理：防止某个梯度方向角度受噪声干扰等因素突变
   • 高斯加权：使特征点附件的梯度幅值具有更大的权重，弥补因没有仿射不变性而导致的特征不稳定
4. 统计出数值最高的梯度方向，作为主方向；保留大于主方向80%的方向作为辅方向

梯度方向和梯度幅值的计算：

$$m(x,y)=\sqrt{(L(x+1,y)-L(x-1,y){)}^2+(L(x,y+1)-L(x,y-1){)}^2}$$

$$\theta (x,y)=\arctan \left({\displaystyle \frac{L(x,y+1)-L(x,y-1)}{L(x+1,y)-L(x-1,y)}}\right)$$

高斯加权公式：

$$W_{i,j}=e^{{\displaystyle \frac{-(i^2+j^2)}{2\times (1.5\sigma {)}^2}}}$$

• $i,j$是像素点到关键点的距离

5. 关键点描述

描述符是一组向量，用于描述特征点及其领域点的特征，以便更好地与其他图片匹配

1. 将坐轴标移到关键点方向
2. 将特征点的领域划分为$d\times d$个子区域，每个子区域大小为$m\sigma \times m\sigma$，并划分为8个方向。论文推荐值：$m=3$，$d=4$
3. 每一个子块进行8个方向的直方图统计操作，共$d\times d\times 8=128$个bin，得到一个128维的特征向量（描述符）

6. 特征点匹配

分别对模板图（参考图，reference image）和实时图（观测图，observation image）建立关键点描述子集合。

目标的识别是通过两点集内关键点描述子的比对来完成，一般采用欧氏距离来衡量两个描述子之间的相似度。

可以使用KD树优化搜索过程。

SIFT实现

cv2中已经给出了SIFT算法的实现，可以直接调用。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

import cv2
import numpy as np
import matplotlib.pyplot as plt

# 读取图像
img=cv2.imread('nz.jpg')
cat=cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)

# 创建SIFT对象
sift=cv2.SIFT_create()

# 关键点检测：kp关键点信息包括方向，尺度，位置信息，des是关键点的描述符
kp,des=sift.detectAndCompute(cat,None)

# 在图像上绘制关键点的检测结果
cv2.drawKeypoints(img,kp,img,flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)

# 图像显示
plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
plt.title('SIFT Result')
plt.show()