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组匹配点对:(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变换的计算,敬请期待!
微信公众号: