图像配准系列之“Sift特征点+薄板样条变换+FFD变换”配准方法

上篇文章中我们讲了“Sift+TPS”的配准方法：

我们知道，TPS薄板样条变换（简称TPS变换）与FFD变换均具有一定的局部形变特性。不过对于TPS变换，只要其输入的匹配点对中有一组发生变化，其整体的形变可能会发生很大改变，因此有时候TPS变换并不能很好地矫正局部形变。而对于FFD变换，每个点的形变坐标偏移只与其周围4*4个控制点的控制参数有关，与其它控制点无关，因此FFD变换的局部形变特性比TPS变换更好，然而“FFD+梯度下降”配准方法的输入参数通常很多，也即要优化的参数很多，优化参数时较容易陷入局部极值，且优化参数的过程很耗时：

图像配准系列之基于FFD形变与梯度下降法的图像配准

基于此原因，有研究员提出了层次FFD的方法，我们前面已经讲过：

图像配准系列之层次FFD形变配准

相比来说“Sift+TPS”的配准方法更稳定更快，为了使配准效果更好、更快、更稳定，因此后来又有人提出进一步把“层次FFD”与“Sift+TPS”结合起来，使配准效果更好、更快、更稳定。如下图所示：

下面上代码：

(1) 层次FFD代码

void level_ffd_match(Mat img1, Mat img2, Mat &outffd)
{
  float min = -0.01;
  float max = 0.01;
  //第一层
  int row_block_num = 8;
  int col_block_num = 8;
  Mat grid_points;
  init_bpline_para(img1, row_block_num, col_block_num, grid_points, min, max);
  Mat out;
  bpline_match(img1, img2, out, row_block_num, col_block_num, grid_points);


  //第二层
  row_block_num = 16;
  col_block_num = 16;
  init_bpline_para(img1, row_block_num, col_block_num, grid_points, min, max);
  Mat out1;
  bpline_match(img1, out, out1, row_block_num, col_block_num, grid_points);




  //第三层
  row_block_num = 30;
  col_block_num = 30;
  init_bpline_para(img1, row_block_num, col_block_num, grid_points, min, max);
  Mat out2;
  bpline_match(img1, out1, out2, row_block_num, col_block_num, grid_points);


  out2.copyTo(outffd);
}

(2) “Sift+TPS+层次FFD”代码

void Sift_Tps_test(void)
{
  Mat img1 = imread("lena.jpg", CV_LOAD_IMAGE_GRAYSCALE);
  Mat img2 = imread("lena_ffd.jpg", CV_LOAD_IMAGE_GRAYSCALE);


  imshow("image before", img1);
  imshow("image2 before", img2);


  // SIFT - 检测关键点并在原图中绘制
  vector<KeyPoint> kp1, kp2;
  kp1 = detect_sift_block(img1, 50, 5, 5);
  kp2 = detect_sift_block(img2, 50, 5, 5);
  Mat outimg1, outimg2;
  drawKeypoints(img1, kp1, outimg1);
  drawKeypoints(img2, kp2, outimg2);
  imshow("image1 keypoints", outimg1);
  imshow("image2 keypoints", outimg2);


  // SIFT - 特征向量提取
  Ptr<SiftDescriptorExtractor> extractor = SiftDescriptorExtractor::create();
  Mat descriptor1, descriptor2;
  extractor->compute(img1, kp1, descriptor1);
  extractor->compute(img2, kp2, descriptor2);


  // 两张图像的特征匹配
  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.35*max_dist)
    {
      goodMatches.push_back(matches[i]);
    }
  }
  cout << "The number of good matches:" << goodMatches.size() << endl;


  //坐标转换为float类型
  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);
  }


  vector<uchar> RANSACStatus;//用以标记每一个匹配点的状态，等于0则为外点，等于1则为内点。
  findFundamentalMat(p01, p02, RANSACStatus, CV_FM_RANSAC, 4.5);//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]);       //基准图特征点
      f2_features_ok.push_back(p02[i]);     //流动图特征点
    }
  }


  Mat Tx, Ty;
  Tps_TxTy(f1_features_ok, f2_features_ok, img2.rows, img2.cols, Tx, Ty);
  Mat tps_out;
  remap(img2, tps_out, Tx, Ty, INTER_CUBIC);  //Sift+TPS粗配准结果


  Mat ffd_out;
  level_ffd_match(img1, tps_out, ffd_out);   //层次FFD细配准结果


  imshow("img2-img1", abs(img2-img1));
  imshow("tps_out-img1", abs(tps_out - img1));
  imshow("ffd_out-img1", abs(ffd_out - img1));
  imshow("tps_out", tps_out);
  imshow("ffd_out", ffd_out);
  cv::waitKey(0);
}

运行上述代码，对扭曲的Lena图像进行配准，结果如下。由以下结果可知，细配准图像比粗配准图像的形变矫正效果好多了。因此结合两种方法的配准方法的配准效果更好。