版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。 本文链接：https://blog.csdn.net/VictoriaW/article/details/70161180 插值 假设已知函数y=f(x)y=f(x)在N+1N+1个点x1,x2,⋯,xN+1x1,x2,⋯,xN+1处的函数值y1,y2,⋯,yN+1y1,y2,⋯,yN+1，但函数的表达式f(x)f(x)未知，那么可以通过插值函数p(x)p(x)来逼近未知函数f(x)f(x)，并且p(x)p(x)必须满足 p(xk)=yk,k=1,2,⋯,N+1.(1) (1)p(xk)=yk,k=1,2,⋯,N+1. 常见的插值函数的形式有多项式函数、样条函数。
多项式函数：令p(x)p(x)为NN次多项式函数，于是p(x)p(x)有N+1N+1个参数，而由公式(1)可知这N+1N+1个参数满足N+1N+1个约束条件，所以可以求出p(x)p(x)的表达式。
样条函数：我们知道NN阶多项式函数必然有N−1N−1个极值点，所以得到的插值函数摆动会比较大，这有点像机器学习中的过拟合现象，可以用样条函数来避免这个问题。这里的样条函数其实就是分段函数，表示在相邻点xkxk和xk+1xk+1之间用低阶多项式函数Sk(x)Sk(x)进行插值。分段线性插值和三次样条插值都属于样条插值。
TPS 本文介绍的TPS针对的是插值问题的一种特殊情况，并且TSP插值函数的形式也比较新颖。 考虑这样一个插值问题：自变量xx是2维空间中的一点，函数值yy也是2维空间中的一点，并且都在笛卡尔坐标系下表示。给定NN个自变量xkxk和对应的函数值ykyk，求插值函数 Φ(x)=[Φ1(x)Φ2(x)]， Φ(x)=[Φ1(x)Φ2(x)]， 使得 yk=Φ(xk).(2) (2)yk=Φ(xk). 我们可以认为是求两个插值函数Φ1(x)Φ1(x)和Φ2(x)Φ2(x)。
TPS插值函数形式如下： Φ1(x)=c+aTx+wTs(x)(3) (3)Φ1(x)=c+aTx+wTs(x) 其中cc是标量，向量a∈ℝ2×1a∈R2×1，向量w∈ℝN×1w∈RN×1，函数向量 s(x)=(σ(x−x1),σ(x−x1),⋯,σ(x−xN))T s(x)=(σ(x−x1),σ(x−x1),⋯,σ(x−xN))T σ(x)=||x||22log||x||2. σ(x)=||x||22log||x||2. Φ2(x)Φ2(x)和Φ1(x)Φ1(x)有一样的形式。看到这里可能会产生疑问？插值函数的形式千千万，怎么就选择公式(3)这种形式呢？我们可以把一个插值函数想象成弯曲一个薄钢板，使得它穿过给定点，这样会需要一个弯曲能量： J(Φ)=∑j=12∬ℝ2(∂2Φj∂x2)2+2(∂2Φj∂x∂y)2+(∂2Φj∂y2)2dxdy J(Φ)=∑j=12∬R2(∂2Φj∂x2)2+2(∂2Φj∂x∂y)2+(∂2Φj∂y2)2dxdy 那么可以证明公式(3)是使得弯曲能量最小的插值函数。参考文献[3]中给了证明过程。
TSP插值函数Φ1Φ1有N+3N+3个参数，而条件(2)只给出了NN个约束，我们再添加三个约束： ∑k=1Nwk=∑k=1Nxxkwk=∑k=1Nxykwk=000(4) (4)∑k=1Nwk=0∑k=1Nxkxwk=0∑k=1Nxkywk=0 xxkxkx和xykxky分别表示点xx的xx坐标值和yy坐标值。于是(2)和(4)可以写成 ⎡⎣⎢⎢⎢S1TNXT1N00X00⎤⎦⎥⎥⎥⎡⎣⎢⎢wca⎤⎦⎥⎥=⎡⎣⎢⎢Yx00⎤⎦⎥⎥(5) (5)[S1NX1NT00XT00][wca]=[Yx00]
其中，(S)ij=σ(xi−xj)(S)ij=σ(xi−xj)，1N1N表示值全为1的NN维列向量， X=⎡⎣⎢⎢⎢⎢xx1xx2⋯xxNxy1xy2⋯xyN⎤⎦⎥⎥⎥⎥ X=[x1xx1yx2xx2y⋯⋯xNxxNy] Yx=⎡⎣⎢⎢⎢⎢yx1yx2⋯yxN⎤⎦⎥⎥⎥⎥ Yx=[y1xy2x⋯yNx] 我们可以令 Γ=⎡⎣⎢⎢⎢S1TNXT1N00X00⎤⎦⎥⎥⎥ Γ=[S1NX1NT00XT00] 那么可知当SS是非奇异矩阵时，ΓΓ也是非奇异矩阵，于是参数为： ⎡⎣⎢⎢wca⎤⎦⎥⎥=Γ−1⎡⎣⎢⎢Yx00⎤⎦⎥⎥(6) (6)[wca]=Γ−1[Yx00] 可以把Φ1Φ1和Φ2Φ2的参数通过一个矩阵运算计算出来： ⎡⎣⎢⎢wxcxaxwycyay⎤⎦⎥⎥=Γ−1⎡⎣⎢⎢Yx00Yy00⎤⎦⎥⎥(7) (7)[wxwycxcyaxay]=Γ−1[YxYy0000] 我们把Γ−1Γ−1写成下面的形式： Γ−1=[Γ11Γ21Γ12Γ22] Γ−1=[Γ11Γ12Γ21Γ22] 称矩阵Γ11Γ11为弯曲能量矩阵，其秩为N−3N−3。
Principal Warp Principal Warp是进一步对TPS进行分析的方法。看原论文[1]的介绍看的好艰辛，等后面如果再碰到的时候再总结。感兴趣的读者可以去阅读论文[1]或者书[2]的第12.3节。
TPS对齐两张图片 假设对齐图片1到图片2。给定图片1和图片2中的已知landmark点，我们可以通过TPS得到由图片1的坐标到图片2的坐标的映射关系。遍历图片1中的每个点，我们可以得到它在图片2中应该对应的点，把图片1中每个点上的像素信息移到对应的图片2中去，就可以得到对准图片2之后的图片1。
但是图片1并不是所有的点都在图片2中有对应，比如：如果图片1中的点映射的横坐标和纵坐标都为负值这种情况。我不知道别人是怎么处理的，目前我是直接舍弃这样的点。
参考 [1] F. L. Bookstein. Principal warps: Thin-plate splines and the decomposition of deformations. IEEE Trans. Pattern Anal. Mach. Intell., 11(6):567–585, June 1989. [2] Ian L. Dryden and Kanti V. Mardia. [Statistical Shape Analysis: With Applications in R].（需要这本书电子版的读者请私信我） [3] Kent, J. T. and Mardia, K. V. (1994a). The link between kriging and thin-plate splines. In: Probability, Statistics and Optimization: a Tribute to Peter Whittle (ed. F. P. Kelly), pp 325–339. John Wiley & Sons, Ltd, Chichester. page 282, 287, 311 ———————————————— 版权声明：本文为博主「Vic时代」的原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接及本声明。 原文链接：https://blog.csdn.net/victoriaw/article/details/70161180
