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

大家好，欢迎来到IT知识分享网。

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组匹配点对：(P1(x1,y1),P1’(x1’,y1’))、P2((x2,y2),P2’(x2’,y2’))、…、(Pn(xn,yn), Pn’(xn’,yn’))。那么使用TPS变换计算图A与图B的坐标对应关系的过程如下。

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

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

上式中只要求出a1，a2，a3和wi(1≤i≤n)，就可以确定f(x,y)，其中U为基函数：

记矩阵K、L、Y为：

由LW=Y解得W矩阵：

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

讲完计算原理，下面我们上基于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变换的计算，敬请期待！