TPS薄板样条变换计算原理及C++实现

TPS薄板样条变换属于一种非刚性形变，该形变算法的输入为两张图像中多组相同部位的匹配点对，输出为两张图像的相同部位的坐标映射，比如图A的点(x1,y1)对应图B的点(x1',y1')，图A的点(x2,y2)对应图B的点(x2',y2')，图A的点(x3,y3)对应图B的点(x3',y3')等等。举个简单例子，假设图B相对于图A的位置，相当于把图A从原来的位置整体向右平移一段距离d，那么经过TPS变换之后得到的坐标对应为：A的每一个像素点向右平移距离d，则是B上对应的点。当然，A与B之间往往还有更复杂的形变，这种情况下坐标对应关系则不是整体一样了，A上每个像素点对应B上每个像素点的情况各有不同，这正体现了TPS变换的非刚性。

TPS变换通常用于图像配准，由于其输入为多组匹配点对，TPS变换之前需要先获取两张图像的多组匹配点对，常见的获取匹配点对算法为Shift、Surf、Orb，以及光流跟踪等。本文主要讲TPS变换的原理，暂时把获取匹配点对的算法放在一边。我们假设已经获取到两张图像的n组匹配点对：(P₁(x₁,y₁),P₁’(x₁’,y₁’))、P₂((x₂,y₂),P₂’(x₂’,y₂’))、…、(P_n(x_n,y_n), P_n’(x_n’,y_n’))。那么使用TPS变换计算图A与图B的坐标对应关系的过程如下。

TPS形变的目标是求解一个函数f，使得f(Pi)=Pi’ (1≤i≤n)，并且弯曲能量函数最小，同时图像上的其它点也可以通过插值得到很好的校正。可以把形变函数想象成弯折一块薄钢板，使这块钢板穿过给定的n个点，弯曲钢板所需能量可表示为：

可以证明薄板样条的插值函数就是弯曲能量最小的函数：

上式中只要求出a₁，a₂，a₃和w_i(1≤i≤n)，就可以确定f(x,y)，其中U为基函数：

记矩阵K、L、Y为：

由LW=Y解得W矩阵：

从而有A的任意坐标(x_i,y_i)到B的任意坐标(x_i’,y_i’)的映射：

讲完计算原理，下面我们上基于C++与Opencv的代码：

基函数的计算：

//TPS基函数，r为两点距离的平方
double tps_basis(double r)
{
  if(r == 0.0)
  {
    return 0.0;
  }
  else
  {
    return (r*log(r));
  }
}

K矩阵的计算：

//n*n矩阵
void cal_K(vector<Point2f> P, Mat &K)
{
  if(K.empty())
  {
    K.create(P.size(), P.size(), CV_32FC1);
  }
  else if(K.rows != P.size() || K.cols != P.size())
  {
    resize(K, K, Size(P.size(), P.size()));
  }
 
 
  for(int i = 0; i < P.size(); i++)
  {
    for(int j = i; j < P.size(); j++)
    {
      Point2f diff_p = P[i] - P[j];
      double diff = diff_p.x*diff_p.x + diff_p.y*diff_p.y;
      float U = (float)tps_basis(diff);
      K.ptr<float>(i)[j] = U;
      K.ptr<float>(j)[i] = U;
    }
  }
}

L矩阵的计算：

//(n+3)*(n+3)矩阵
void cal_L(vector<Point2f> points1, Mat &L)
{
  int N = points1.size();
  
  Mat L_tmp =  Mat::zeros(Size(N+3, N+3), CV_32FC1);
  
  Mat P = Mat::ones(Size(3, N), CV_32FC1);
  for(int i = 0; i < N; i++)
  {
    P.ptr<float>(i)[1] = points1[i].x;
    P.ptr<float>(i)[2] = points1[i].y;
  }
 
 
  Mat P_T;
  transpose(P, P_T);
 
 
  Mat K;
  cal_K(points1, K);
 
 
  K.copyTo(L_tmp(Rect(0, 0, N, N)));
  P.copyTo(L_tmp(Rect(N, 0, 3, N)));
  P_T.copyTo(L_tmp(Rect(0, N, N, 3)));
 
 
  L = L_tmp.clone();
}

W矩阵的计算：

//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
}

坐标映射的计算，计算图A与图B中每一个点的坐标映射关系：

Point2f Tps_Transformation(vector<Point2f> points1, Point2f point, Mat W)
{
 
 
    Point2f out;
 
 
    float a1_x = W.at<float>(W.rows-3, 0);
    float ax_x = W.at<float>(W.rows-2, 0);
    float ay_x = W.at<float>(W.rows-1, 0);
 
 
    float a1_y = W.at<float>(W.rows-3, 1);
    float ax_y = W.at<float>(W.rows-2, 1);
    float ay_y = W.at<float>(W.rows-1, 1);
 
 
    float affine_x = a1_x+ax_x*point.x + ay_x*point.y;
    float affine_y = a1_y+ax_y*point.x + ay_y*point.y;
 
 
    float nonrigid_x = 0;
    float nonrigid_y = 0;
 
 
    for (int j = 0; j < points1.size(); j++)
    {
         Point2f diff_p = points1[j] - point;
         double diff = diff_p.x*diff_p.x + diff_p.y*diff_p.y;
         double tps_diff = tps_basis(diff);
         nonrigid_x += W.at<float>(j, 0)*tps_diff;
         nonrigid_y += W.at<float>(j, 1)*tps_diff;
    }
 
 
    out.x = affine_x + nonrigid_x;
    out.y = affine_y + nonrigid_y;
    
    return out;   //返回(x, y)对应的点坐标(x', y')
}

由于每个点的坐标分为x坐标与y坐标，因此把获取到的点应该拆分为坐标与y坐标，所以最后得到Tx与Ty两个Mat矩阵，分别对应所有像素点的x坐标映射与y坐标映射。

//获取图A与图B的每个点的映射坐标
void Tps_TxTy(vector<Point2f> points1, vector<Point2f> points2, int row, int col, Mat &Tx, Mat &Ty)
{
 
 
    Mat mapX(row, col, CV_32FC1);
    Mat mapY(row, col, CV_32FC1);
 
 
 
 
    Mat W;
    cal_W(points1, points2, W);
 
 
    for (int i = 0; i < row; i++)
    {
        for (int j = 0; j < col; j++)
        {
            Point2f pt = Tps_Transformation(points1, Point2f(float(j), float(i)), W);
            mapX.at<float>(i, j) = pt.x;
            mapY.at<float>(i, j) = pt.y;
        }
    }
 
 
  mapX.copyTo(Tx);
  mapY.copyTo(Ty);
}

以上代码运行起来是非常慢的，因为其需要计算对数，每一个点的坐标映射大约需要计算n次（n组匹配点）对数与惩罚，耗时点主要在这里。在之后的文章中，我们再讲解使用CUDA来加速TPS变换的计算，敬请期待！

微信公众号：