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

在前文中，我们已经有好几篇讲图像配准的文章：

在以上第1、第2篇文章中，我们讲了图像配准的概念、Sift特征点的基本原理，以及基于“Sift特征点+仿射变换”的图像配准原理和C++实现。我们知道，仿射变换模型具有整体的特性，对局部形变的矫正效果并不理想。相比来说FFD变换与TPS薄板样条变换，则对局部形变具有更强的矫正能力。基于FFD变换的配准方法我们在前文已经讲过，本文我们主要讲基于“Sift特征点+TPS薄板样条变换”的配准方法，并使用C++与Opencv来将之实现。

1. Sift算法原理

在上述第1、2篇文章中已讲，此处不再重复。

2. TPS薄板样条变换原理

在上述第8篇文章中已讲，此处也不再重复。

3. 配准步骤

配准步骤的大致步骤如下图所示：

其中特征向量的匹配与剔除，我们在上方超链接的第2篇文章已讲过。

像素重采样时，由于坐标是浮点型的，会使用到插值算法，我们在之前的文章也详细说明：

常见图像插值算法的原理与C++实现

双三次插值算法的C++实现与SSE指令优化

4. 代码实现

好了，原理我们基本都讲过，下面直接上代码。

(1) TPS的代码

可参考开头超链接的第8篇文章，不过这里需要修正一下矩阵W的计算，调用Opencv的solve函数解矩阵方程LW=Y时，如果参数是DECOMP_LU时会出错，但是如果改为DECOMP_EIG就没问题，这里应该涉及到矩阵可逆以及求解算法的问题，有空再继续深究。

//W = {w0, w1, w2, ..., a1, ax, ay}
void cal_W(vector<Point2f> points1, vector<Point2f> points2, Mat &W)
{
  int N = points1.size();
  Mat L;
  cal_L(points1, L);


  Mat Y = Mat::zeros(Size(2, N+3), CV_32FC1);  //row:N+3, col:2
  for(int i = 0; i < N; i++)
  {
    Y.ptr<float>(i)[0] = points2[i].x;
    Y.ptr<float>(i)[1] = points2[i].y;
  }


  //solve(L, Y, W, DECOMP_LU);   //LW = Y,求W
  solve(L, Y, W, DECOMP_EIG);
  //W = L.inv()*Y;
}

(2) 检测特征点代码

这里有个问题需要说明一下，如果直接对整张图像进行Sift特征点检测，得到的特征点分布不均匀，然而点分布越均匀，TPS变换效果越好，因此这样直接检测整张图的特征点，最后配准效果并不理想。所以我们将图像分为多个块，分别检测各个块的特征点，然后拼接成整张图像的特征点，这样一来，特征点的分布就比较均匀了。

代码如下：

vector<KeyPoint> detect_sift_block(Mat img, int num, int b_row, int b_col)
{
  const int block_rsize = img.rows / b_row;
  const int block_csize = img.cols / b_col;


  vector<KeyPoint> kp;
  vector<KeyPoint> kp_tmp;
  Ptr<SiftFeatureDetector> siftdtc = SiftFeatureDetector::create(num, 3, 0.04, 10.0, 1.6);
  
  for (int i = 0; i < b_row; i++)   //分块检测，确保均匀分布
  {
    for (int j = 0; j < b_col; j++)
    {
      int xx = j*block_csize;
      int yy = i*block_rsize;
      Mat tmp;
      img(Rect(xx, yy, block_csize, block_rsize)).copyTo(tmp);
      siftdtc->detect(tmp, kp_tmp);
      
      //将小块的点坐标转换成整张图像的点坐标
      for (int k = 0; k < kp_tmp.size(); k++)
      {
        kp_tmp[k].pt.x += xx;
        kp_tmp[k].pt.y += yy;
      }
      //将多个小块的特征点拼接成整张图像的特征点
      kp.insert(kp.end(), kp_tmp.begin(), kp_tmp.end());
    }
  }


  return kp;
}

(3) 配准代码

使用具有局部形变的Lena图像进行配准测试，分别使用TPS变换与仿射变换作为坐标变换模型，并比较结果。

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]);     //流动图特征点
    }
  }


#if 1   //这里使用TPS变换
  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);


  imshow("img2-img1", abs(img2-img1));
  imshow("tps_out-img1", abs(tps_out - img1));
  imshow("tps_out", tps_out);
#else  //这里使用仿射变换
  Mat T = estimateRigidTransform(f2_features_ok, f1_features_ok, true);   //false表示刚性变换


  Mat affine_out;
  warpAffine(img2, affine_out, T, img2.size(), INTER_CUBIC);


  imshow("img2-img1", abs(img2 - img1));
  imshow("affine_out-img1", abs(affine_out - img1));
  imshow("affine_out", affine_out);
#endif


  cv::waitKey(0);
}

5. 配准结果

(1) 首先是使用“Sift+仿射变换”的结果：