特征匹配：GMS(2017): Grid-based Motion Statistics for Fast, Ultra-robust Feature Correspondence

Reference:

1. GMS: Grid-based Motion Statistics for Fast, Ultra-robust Feature Correspondence
2. GMS特征匹配 原文和代码详细解读GMS: Grid-based Motion Statistics for Fast, Ultra-robust Feature Correspondence
3. 特征匹配算法GMS2020(Grid-based Motion Statistics)理论与实践
4. GMS: Fast and Robust Feature Matcher (CVPR 17 & IJCV 20)

该方法现已集成进 OpenCV 库 xfeatures2d.matchGMS。

1. 简介

将平滑约束纳入特征匹配以实现超鲁棒匹配。然而，这样的公式既复杂又缓慢，不适合视频应用。本文提出了基于网格的运动统计(GMS - Grid-based Motion Statistics)，这是一种将运动平滑度封装为区域中一定数量的匹配的统计可能性的简单方法。GMS 能够将高匹配数转换成高匹配质量。这提供了实时、超鲁棒对应系统。对纹理低、模糊和宽基线的视频的评估表明，GMS始终优于其他实时匹配器，并且可以与更复杂、更慢的技术达到相似的效果。

与 SIFT 相比效果如下：
本文主要贡献如下：

1. 将运动平滑约束转换为拒绝错误匹配的统计度量；
2. 提出一个高效的基于网格的分数估计器，可以整合到实时特征匹配器中；
3. 将 GMS 方法与传统的 SIFT、SURF、LITF 特征进行对比，表明 GMS 方法在时间与准确率方面都要优于其他算法。

1.1 效果评估

视频中截取的效果对比如下：
GMS 算法与其他算法在不同数据集上的性能与耗时对比如下。其中横轴表示时间，纵轴表示性能：

可以看到 GMS 确实实现了高精度、低耗时的特征匹配。

2. 方法

2.1 符号

文章主要用到的符号定义如下：

符号	解释
$x_i$	第 $i$ 个特征匹配
$f_a$	区域 $a$ 中的 $n$ 个支持特征中的一个
$S_i$	特征匹配 $x_i$ 的邻域评分
{ p t , p f } \{p_t, p_f\} {pt​,pf​}	分别代表 { p ( f a b ∣ T a b ) , p ( f a b ∣ F a b ) } \{p(f_a^b|T^{ab}), p(f_a^b|F^{ab})\} {p(fab​∣Tab),p(fab​∣Fab)}
$T^{ab}$	区域 { a , b } \{a, b\} {a,b}看到相同的位置
$F^{ab}$	区域 { a , b } \{a, b\} {a,b}看到不同的位置
$f_a^t$	若 $f_a$ 匹配正确，则 $p\left(f_a^t\right)=\mathrm{t}$
$f_a^f$	若 $f_a$ 匹配错误，则 $p\left(f_a^f\right)=1-\mathrm{t}$
$f_a^b$	$f_a$ 的最近邻匹配特征在区域 $b$ 内(即$x_i$的邻居的最近邻匹配如果在区域$b$内，该概率增加)
$\overline{f_a^b}$	$f_a$ 的最近邻匹配特征不在区域 $b$ 内
$m$	在区域 $b$ 内的特征个数
$n$	在区域 $a$ 内的特征个数

1.2 论文核心思路

文章核心就一句话：正确匹配周围会有较多的匹配去支持它，而错误的匹配周围支持它的匹配很少。

• 假设一： 如果运动是平滑的，相邻像素和特征应该一起移动。因此就有了以下假设：如果两个特征点匹配正确，那么运动平滑度导致真实匹配周围的(小)邻域可以看做对应同一个3D位置(原文中写的就是location，但我感觉它的意思还是区域)；而错误匹配的邻域对应不同的3D位置。这里邻域被定义为围绕下图所示的相应图像特征的一对区域 { a , b } \{a,b\} {a,b}：
  
  从图中可以看出，对于一个匹配 $x_i$，在 $I_a$ 和 $I_b$ 两幅图上的区域分别记为 { a , b } \{a, b\} {a,b}。若 $x_i$ 匹配正确，则 $a$ 和 $b$ 对应同一个 3D 区域，所以 $a$ 和 $b$ 内匹配的点对数会很多；同理，若 $x_j$ 匹配错误，则对应的小邻域内匹配的点对数就很少。引入一个分数估计方法 $S_i$，表示匹配点 $x_i$ 小邻域内匹配的点对数，上图中若 $x_i$ 邻域内的匹配点对数为 $2$(黄框)，则 $S_i=2$，$x_j$ 邻域内的匹配点个数为 $0$，则 $S_j=0$(红框)。

• 假设二：如果给定 $f_a$ 匹配错误，它的最近邻正确匹配可以位于 $I_b$ 所有特征点中的任何一个，即在整幅图片所有的$M$个特征点内的概率相同。之所以有这个假设，是因为一般来说，特征的错误最近邻是因为没有先验原因来支持任何图像区域。
  假设在 $I_a$ 和 $I_b$ 两幅图上分别有 $N$ 和 $M$ 个特征点，则有(在 $f_a$ 匹配错误的条件下，$f_a$ 的最近邻匹配特征在区域 $b$ 内的概率)：
  $$p\left(f_a^b\mid f_a^f\right)=\beta m/M$$其中 $m$ 表示在区域 $b$ 内的特征个数，$\beta$ 是添加的因子，用它来适应如一行窗户等重复的结构而违反该假设(一堆一样的窗户，很容易匹配错)。

  给定 { a , b } \{a,b\} {a,b} 能看到相同的位置，特征 $f_a$ 的最近邻居在区域 $b$ 内，则：
  $$\begin{array}{rl}{p}_t=p\left(f_a^b\mid T^{ab}\right)& =p\left(f_a^t\mid T^{ab}\right)+p\left(f_a^f,f_a^b\mid T^{ab}\right)\\ & =p\left(f_a^t\mid T^{ab}\right)+p\left(f_a^f\mid T^{ab}\right)p\left(f_a^b\mid f_a^f,T^{ab}\right)\\ & =p\left(f_a^t\right)+p\left(f_a^f\right)p\left(f_a^b\mid f_a^f\right)\\ & =t+(1-t)\beta m/M\end{array}$$下图表示了 $f_a^b$ 仅当 $f_a$ 正确匹配(下图绿色)或错误匹配但巧合的落在区域 $b$ 中(下图中间的柠檬黄色)时会发生。这给出了上面等式的第一行，即 $p\left(f_a^t\mid T^{ab}\right)+p\left(f_a^f,f_a^b\mid T^{ab}\right)$，第二行是贝叶斯规则。由于特征是预匹配的，$p(f_a^t)$、$p(f_a^f)$ 与 $T^{ab}$ 无关。由于假设2，$p\left(f_a^b\mid f_a^f\right)$ 也与 $T^{ab}$ 无关，最终得到上式最终的表达。

Figure 3(i):

• 绿：理论上应该是匹配到 b 区域的对应特征点，也匹配正确了的邻居概率；
• 柠檬黄：理论上应该是匹配到 b 区域的对应特征点，匹配错了特征点但依旧在 b 区域内；
• 橙：理论上应该是匹配到 b 区域的对应特征点，匹配错了且特征点不在 b 区域内。

Figure 3(ii):

• 绿：理论上应该匹配到不在 b 区域的特征点，这时正确匹配上了；
• 柠檬黄：理论上应该匹配到不在 b 区域的特征点，这时匹配错点了，而很不凑巧误匹配到了区域 b 中的特征点；
• 橙：理论上应该匹配到不在 b 区域的特征点，这时匹配错点了，但是它不在区域 b 中。

令 $p_f=p\left(f_a^b\mid F^{ab}\right)$，与上面公式类似，可以得到：
$$\begin{array}{rl}{p}_f=p\left(f_a^b\mid F^{ab}\right)& =p\left(f_a^f,f_a^b\mid F^{ab}\right)\\ & =p\left(f_a^f\mid F^{ab}\right)p\left(f_a^b\mid f_a^f,F^{ab}\right)\\ & =p\left(f_a^f\right)p\left(f_a^b\mid f_a^f\right)\\ & =\beta (1-t)(m/M)\end{array}$$由于每个特征的匹配是独立的，使用 假设1 和上面两个等式，我们可以近似 $S_i$ 分布------一对二项分布来表示 $x_i$ 邻域内的匹配数：
$$S_i\sim \{\begin{array}{ll}B\left(n,p_t\right),& \text{ if }{x}_i\text{ is true }\\ B\left(n,p_f\right),& \text{ if }{x}_i\text{ is false }\end{array}$$这里的要点是真匹配和误匹配具有的领域分数 $S$ 有着非常不同的分布，如下图所示，所以这里的分数 $S$ 可以称为区分真假匹配的有用指标。

运动在大范围内通常是平滑的。然而，假设1 需要足够小的邻域。在超大邻域中，真正的匹配邻域将包含一些错误匹配区域，反之亦然。这减少了真假分数分布的可分离性。因此，概括 假设1，我们有：

• 假设三：如果运动在区域上平滑，则真正的匹配允许预测查看相同 3D 位置的多个小区域对。在错误匹配上使用相同的预测函数将导致几何上不同的3D位置。也就是说，可以将一个较大的区域分成 $K$($K$不相交) 个很小的区域，再次用以上公式，二项分布就变成了：
  $$S_i\sim \{\begin{array}{ll}B\left(Kn,p_t\right),& \text{ if }{x}_i\text{ is true }\\ B\left(Kn,p_f\right),& \text{ if }{x}_i\text{ is false }\end{array}$$$S_i$ 分布的均值和方差分别为：
  $\left\{m_t=Kn{p}_t,s_t=\sqrt{Knt\left(1-p_t\right)}\text{ if }{x}_i\text{ is true}\right\}$
  $\left\{m_f=Kn{p}_f,s_f=\sqrt{Kn{p}_f\left(1-p_f\right)}\text{ if }{x}_i\text{ is false}\right\}$

  在统计事件中，如果一个事件偏离其中值很多个标准差我们会认为此事件不会发生，我们因此会希望以下公式越大越好，$P$ 越大就越具有区分度：
  $$P=\frac{m_t-m_f}{s_t+s_f}=\frac{Kn{p}_t-Kn{p}_f}{\sqrt{Kn{p}_t\left(1-p_t\right)}+\sqrt{Kn{p}_f\left(1-p_f\right)}}.$$

---

增大 $P$ 有两种方式：增大 $K$ 或者增大 $n$。传统的特征匹配方法侧重于增加描述符不变性，由下图可见：

在上面一行中，大变形的 A-SIFT 具有最佳的宽基线性能。而 GMS 使用标准ORB通过简单地增加特征点数量(括号中数量)，实现了与上面相同的效果。这使得ORB成为一个可行的宽基线特性。

可配置项： { m , n , β , t , K , G } \{m,n,\beta,t,K,G\} {m,n,β,t,K,G}

2.1 加速匹配原理

2.1.1 通过网格单元进行有效分数计算

将图像划分为 $G=20\times 20$ 的不重叠单元格，每个单元格对只计算一次(也就是说每个块计算一次评分？而不是每个点单独算评分)，这样计算每个特征邻域的时间复杂度从 $O(N)$ 砍到了 $O(1)$。在实践中，许多特征位于网格边缘，因此多重复三次网格计算，分别将 $x$ 和 $y$ 平移半个单元格。

若遇到更加复杂纹理多的图片可以增加网格个数。

2.1.2 哪些邻域（网格单元）组合在一起

划分出的大网格是不具备足够小这一假设的，所以每个大网格还要继续划分，令 $K=9$，即每个大网格还要继续划分成 $9$ 个小网格，所以对单元格对 { i , j } \{i, j\} {i,j} 的评分 $S_{ij}$ 的方程为：
$${\mathcal{S}}_{ij}=\sum _{k=1}^{K=9}\left|{\mathcal{X}}_{i^k{j}^k}\right|$$其中 $\left|{\mathcal{X}}_{i^k{j}^k}\right|$ 是上图单元格 { i k , j k } \{i^k, j^k\} {ik,jk} 间的匹配数。

这种方案有效但限制了平面内旋转不变性，虽然依旧比较有竞争力，但是明显在较低旋转下效果更好(注意在旋转90～180的情况下，效果就很差了)。

3. 代码详细讲解

1. gms_matcher gms(kp1, img1.size(), kp2, img2.size(), good_matches); 构建gms_matcher对象，在默认构造器中：
   ① NormalizePoints(vkp1, size1, mvP1); 归一化点的坐标到(0,1)区间；
   ② ConvertMatches(vDMatches, mvMatches); 把opencv的matches转换成vector形式的matches；
   ③ InitalizeNiehbors(mGridNeighborLeft, mGridSizeLeft); mGridNeighborLeft 为 $400\ast 9$ 的Mat，$400$ 对应 $20\ast 20$ 块，$9$ 对应每一块的邻块的序号。
2. gms.GetInlierMask(vbInliers, false, false); 开始做GMS匹配筛选，第一个false表示不做多尺度(代码定义了5种尺度mScaleRatios)缩放筛选，第二个false表示不做多旋转方向上(代码定义了8种旋转mRotationPatterns)的筛选，筛选后的结果保存在vbInliers中.
   ① void SetScale(int Scale) 对不同scale，缩放图片大小，再InitalizeNiehbors(mGridNeighborLeft, mGridSizeLeft); 这个函数参照上面1中的③。
3. int gms_matcher::run(int RotationType) 正式进入筛选，参数RotationType为1~8的旋转编号。
    ① AssignMatchPairs(int GridType) 该函数主要是为了统计两个参数-----Mat mMotionStatistics 和 vector<int> mNumberPointsInPerCellLeft。
     对于前者，mMotionStatistics.at<int>(lgidx, rgidx) 的值表示，在区块lgidx和rgidx中的匹配对的数目；
     对于后者，mNumberPointsInPerCellLeft[lgidx]的值表示区块lgidx和右图所有的匹配对的数目，不一定在rgidx中。
    ② VerifyCellPairs(int RotationType) 核心函数，判断匹配的区块是否是inlier块。
     1) 遍历左图每一区块，对每一区块 $i$，遍历右图区块找到与 $i$ 匹配对最多的块 $j$；
     2) 同时遍历 $i$ 和 $j$ 的 3 * 3 邻块，统计每对邻块的匹配点数；
     3) 计算阈值，${\tau }_i=\alpha \sqrt{n_i}$，滤除小于阈值的点，其中$\alpha$ 经验值为 6，$n_i$ 表示第 $i$ 区块加上3 * 3邻块的平均每块的匹配数量。

4. 代码使用

std::vector<cv::DMatch> good_matches;
// 汉明距离匹配方式
cv::BFMatcher matcher(cv::NORM_HAMMING);
matcher.match(kf1->mDescriptors, kf2->mDescriptors, good_matches);

std::vector<bool> inlinerMask;
GMSMatcher test_gms(kf1->mvKeys, kf1->mImGray.size(), kf2->mvKeys, kf2->mImGray.size(), good_matches);
test_gms.GetInlierMask(inlinerMask, false, false);

std::vector<cv::DMatch> final_matches;
for (size_t index = 0; index < inlinerMask.size(); index++) {
  if (inlinerMask[index]) {
    final_matches.push_back(good_matches[index]);
    matches.emplace_back(good_matches[index].queryIdx, good_matches[index].trainIdx);
  }
}
cv::Mat img_matches_gms;
drawMatches( kf1->mImGrays[cam_id], kf1->mvKeys[cam_id], kf2->mImGrays[cam_id], kf2->mvKeys[cam_id], final_matches, img_matches_gms );
imshow("Matches_gms", img_matches_gms );