图像配准系列之Sift特征点检测与匹配(2)

在上一篇文章中，我们讲到了梯度直方图的抛物线插值，至此，我们得到了每个极值点的位置（x，y坐标）、尺度（σ）、方向（θ）。接下来则是利用这些信息来构建极值点的128维向量描述符，以及描述符的匹配、误匹配的剔除、变换矩阵的计算、像素重映射等。

https://blog.csdn.net/shandianfengfan/article/details/113156703?spm=1001.2014.3001.5501

1. 根据位置、尺度、方向信息构建特征点描述符

(1) 确定计算区域

尺度决定了极值点对应于高斯金字塔图像中的哪个大层和哪个小层的图像Lσ，极值点的描述符就是描述在Lσ上该点及其邻域的梯度信息。如下图，思路是取以极值点为中心的d*d(通常取d=4)个小正方形区域来计算描述符，而每个小正方形的边长为3σ，故计算区域是边长为3σ*d的正方形。

由于生成梯度直方图需要进行三线性插值操作，在以上思路的基础上，计算区域变为3σ*(d+1)的正方形。

如上图，为了使描述子具有旋转不变性，需要把计算描述符的正方形区域旋转θ角度，这样一来即使匹配的特征点在图像中角度不一样，但它们各自计算描述子的区域都已经旋转到各自的方向上面，所以计算所得描述子体现的方向是一致的。正方形区域的外接圆半径为：

为了取得旋转之后正方形区域的所有像素点，实际计算时我们把计算区域扩大为边长是D的外接正方形（上图中外接圆的外接正方形），D的计算如下式：

以极值点为中心、以D为边长的外接正方形，就是该极值点的描述符计算区域，但实际对描述符有贡献的区域则是其内部经过旋转的正方形区域。假设旋转之前内部正方形区域的任意点坐标为(x,y)，旋转之后对应点的坐标为(x',y')，那么(x',y')可按下式计算：

(2) 计算描述符

经过以上步骤，我们得到了计算描述符的区域，接下来还需分别统计旋转之后的4*4个小正方形区域的梯度直方图。按照本人的理解，把0~360°均分成8个角度区间，对于每个小正方形区域中的每个点，其梯度方向属于哪个区间，则把其梯度幅值累加到该区间的直方图幅值，共有4*4个小正方形，每个小正方形有8个区间，每个区间对应1个直方图幅值，，因此把4*4*8=128个梯度直方图幅值构成的128维向量，则是最终的特征点：极值点+描述符。

然后实际上，为了使特征点更加稳定，还增加了以下两个步骤：

I. 高斯加权。以σ=0.5d作为标准差，计算小正方形中每个点的权重，并把该点的梯度幅值乘以高斯权重，作为最终的梯度幅值。

II. 三线性插值。这个本人还没能理解透彻，到底怎么个插值法，改天再深入研究三线性插值原理以及在Sift算法中的应用。

III. 得到每个极值点的128维特征向量之后，为了消除光照亮度差异的影响，还需要对向量进行归一化。假设128维向量为h1、h2、h3、...、h128，按下式对其归一化，即可得到最终的描述符H'。

2. 特征向量的匹配

假设图A检测到m个特征点，图B检测到n个特征点，那么计算出图A与图B的两两特征点的相似度，对于图A的每个特征点P，在图B中与其最相似的特征点是P'，则认为P与P'为图A与图B中互相匹配的特征点。向量相似度的衡量通常使用欧式距离：

初步得到两图中多组互相匹配的点对之后，往往还存在不少误匹配点对，这时还需要将这些误匹配的点对剔除：

(1) 假设所有匹配点对中，距离最大值为max，判断每组点对的距离是否满足下式，如果不满足则将其剔除，其中α取一个小于1的数，比如0.25。

(2) 使用RANSAC算法进一步剔除误匹配点对。

RANSAC算法的一般步骤如下：

a. 随机从数据集中抽取计算数学模型参数所需的最少样本（抽取的样本不能共线），计算模型参数，记所得参数确定的模型为M；

b. 由第1步得到模型M之后，使用该模型计算所有数据的投影误差，并将误差小于阈值的数据包含到内点集I中；

c. 假如I中包含的数据量大于历史最优内点集I_best，那么使用I包含的数据完全替换掉I_best中的数据，并更新最大迭代次数k，k随I_best的变化而改变。其计算公式如下式所示，其中w为I_best的内点比例，p为置信度，通常取0.995，m为计算模型参数所需的最少样本数；

d. 当迭代次数大于k则退出迭代，否则将迭代次数加1并重复以上步骤。

使用RANSAC算法剔除异常匹配点对时，通常选择3×3的单应矩阵作为数学模型，并令矩阵的最后一个参数为固定值1来归一化矩阵，因此实际上单应矩阵模型只需要确定8个未知参数即可。假设随机抽取参考图像上的一个待跟踪点坐标为(xs’,ys’)，浮动图像上的对应跟踪点坐标为(xs,ys)，那么有如下关系（其中s为尺度参数）：

由上式可得：

1组匹配点对可以成立上式的两2个方程，因此需要随机抽取4组不共线的匹配点对来确定单应矩阵中的8个未知参数。单应矩阵的参数确定之后，再使用单应矩阵计算所有匹配点对的投影误差，如下式：

根据上式中的投影误差f是否小于设定阈值来判断匹配点对是否为合理匹配，如果小于则认为合理，否则认为是异常匹配。按上述四个步骤重复迭代计算单应矩阵的参数，直到迭代次数超过最大迭代次数为止，从而找到并剔除异常匹配点对。

3. 坐标变换的计算

坐标变换模型通常有投影变换、仿射变换、薄板样条变换、基于B样条的自由形变等，这里我们只讲最常使用的变换模型之一：仿射变换。

仿射变换矩阵为3×3矩阵，但由于其最后一行为固定的0,0,1，因此通常将其简化为2×3矩阵。假设图A与图B的一组匹配点对坐标分别为(x, y)和(x’, y’)，则有以下计算关系：

由式上式得到：

图A与图B的一系列匹配点对分别记为(P1(x1,y1),P1’(x1’,y1’))、P2((x2,y2),P2’(x2’,y2’))、…、(Pn(xn,yn), Pn’(xn’,yn’))。可以通过解超定方程组来求得满足最小二乘法的仿射变换参数。

我们记：

只需解方程AX=B得到矩阵X即可确定仿射变换参数。由AX=B得到：

最后通过计算式上式的X，即可得到仿射变换的6个参数。

4. 像素重采样

像素重采样在之前的文章我们有讲过，感兴趣的读者可参阅以下博文：

https://blog.csdn.net/shandianfengfan/article/details/111570598?spm=1001.2014.3001.5501

好了，基于Sift算法的图像配准原理我们就暂时讲到这里，其原理还有待进一步深究。下面我们上代码，讲解Opencv中的Sift算法模块怎么使用。

#include <iostream>
#include <opencv2/opencv.hpp>
#include <opencv2/core/core.hpp>
#include <opencv2/xfeatures2d.hpp>
#include <opencv2/features2d/features2d.hpp>


using namespace cv;
using namespace std;
using std::cout;
using std::endl;
using std::vector;
using cv::Mat;
using cv::xfeatures2d::SiftFeatureDetector;
using cv::xfeatures2d::SiftDescriptorExtractor;


void Sift_test(void)
{
  // 从文件中读入图像
  Mat img1 = imread("img1.jpg", CV_LOAD_IMAGE_GRAYSCALE);
  Mat img2 = imread("img2.jpg", CV_LOAD_IMAGE_GRAYSCALE);
  
  imshow("image before", img1);
  imshow("image2 before", img2);


  // SIFT - 检测关键点并在原图中绘制
  int kp_number{ 200 };   //设定最多检测200个特征点
  vector<KeyPoint> kp1, kp2;
  //创建特征点检测对象
  Ptr<SiftFeatureDetector> siftdtc = SiftFeatureDetector::create(kp_number);
  siftdtc->detect(img1, kp1);  //检测图1的极值点
  Mat outimg1;
  drawKeypoints(img1, kp1, outimg1);  //画出图1中检测到的极值点
  imshow("image1 keypoints", outimg1);
  
  siftdtc->detect(img2, kp2);   //检测图2的极值点
  Mat outimg2;
  drawKeypoints(img2, kp2, outimg2);  //画出图2中检测到的极值点
  imshow("image2 keypoints", outimg2);


  // SIFT - 特征向量提取
  Ptr<SiftDescriptorExtractor> extractor = SiftDescriptorExtractor::create();
  Mat descriptor1, descriptor2;
  extractor->compute(img1, kp1, descriptor1);  //计算图1的极值点地描述符
  extractor->compute(img2, kp2, descriptor2);  //计算图2的极值点地描述符


  //创建特征点匹配对象
  Ptr<DescriptorMatcher> matcher = DescriptorMatcher::create("BruteForce");
  vector<DMatch> matches;
  Mat img_matches;
  matcher->match(descriptor1, descriptor2, matches);  // 两张图像的特征匹配
  
  //计算匹配结果中距离最大和距离最小值
  double min_dist = matches[0].distance, max_dist = matches[0].distance;
  for (int m = 0; m < matches.size(); m++)
  {
    if (matches[m].distance < min_dist)
    {
      min_dist = matches[m].distance;
    }
    if (matches[m].distance > max_dist)
    {
      max_dist = matches[m].distance;
    }
  }
  cout << "min dist=" << min_dist << endl;
  cout << "max dist=" << max_dist << endl;


  
  vector<DMatch> goodMatches;
  for (int i = 0; i < matches.size(); i++)   //筛选出较好的匹配点对
  {
    if (matches[i].distance < 0.25*max_dist)  //小于最大距离一定比例的点对则认为合格，否则剔除
    {
      goodMatches.push_back(matches[i]);
    }
  }
  cout << "The number of good matches:" << goodMatches.size() << endl;


  Mat good_img_matches;
  drawMatches(img1, kp1, img2, kp2, goodMatches, good_img_matches);
  imshow("good_img_matches", good_img_matches);




  //将以上初步剔除时满足条件的特征点对筛选出来
  vector <KeyPoint> good_kp1, good_kp2;
  for (int i = 0; i < goodMatches.size(); i++)   
  {
    good_kp1.push_back(kp1[goodMatches[i].queryIdx]);
    good_kp2.push_back(kp2[goodMatches[i].trainIdx]);
  }


  //获取匹配特征点对的坐标值
  vector <Point2f> p01, p02;
  for (int i = 0; i < goodMatches.size(); i++)
  {
    p01.push_back(good_kp1[i].pt);
    p02.push_back(good_kp2[i].pt);
  }


  //RANSAC算法进一步剔除误匹配点对
  vector<uchar> RANSACStatus;//用以标记每一个匹配点的状态，等于0则为外点，等于1则为内点。
  findFundamentalMat(p01, p02, RANSACStatus, CV_FM_RANSAC, 3);//p1 p2必须为float型
  vector<Point2f> f1_features_ok;
  vector<Point2f> f2_features_ok;
  for (int i = 0; i < p01.size(); i++)   //剔除跟踪失败点
  {
    if (RANSACStatus[i])
    {
      f1_features_ok.push_back(p01[i]);       //图1的点
      f2_features_ok.push_back(p02[i]);     //图2的对应点
    }
  }


  //计算仿射变换矩阵
  Mat T = estimateRigidTransform(f2_features_ok, f1_features_ok, true);   //false表示刚性变换


  Mat affine_out;
  //像素重采样
  warpAffine(img2, affine_out, T, img2.size(), INTER_CUBIC);


  imshow("affine_out", affine_out);


  waitKey(0);
}

运行上述代码，对具有平移和旋转的图像进行配准，结果如下：

好了，基于Sift算法的图像配准我们就暂时讲到这里，由于本人水平有限，难免有理解不当的地方，欢迎各位指正。