在上一篇文章中，我们讨论了梯度直方图的抛物线插值。到当前为止，我们已经获得了每个极值点的位置(X，Y坐标)、标度()和方向()。接下来利用这些信息构造极值点的128维向量描述子，以及描述子的匹配、误匹配的消除、变换矩阵的计算、像素重映射等。 　　

　　

　　https://blog . csdn . net/shandian fengfan/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 #include #include #include #include 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 kp1, kp2;  //创建特征点检测对象  Ptr 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 extractor = SiftDescriptorExtractor::create();  Mat descriptor1, descriptor2;  extractor->compute(img1, kp1, descriptor1);  //计算图1的极值点地描述符  extractor->compute(img2, kp2, descriptor2);  //计算图2的极值点地描述符  //创建特征点匹配对象  Ptr matcher = DescriptorMatcher::create("BruteForce");  vector 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 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  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  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 RANSACStatus;//用以标记每一个匹配点的状态，等于0则为外点，等于1则为内点。  findFundamentalMat(p01, p02, RANSACStatus, CV_FM_RANSAC, 3);//p1 p2必须为float型  vector f1_features_ok;  vector 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算法的图像配准我们就暂时讲到这里，由于本人水平有限，难免有理解不当的地方，欢迎各位指正。