【1】matlab程式碼https://www.mathworks.com/matlabcentral/fileexchange/24315-warping-using-thin-plate-splines

【2】

url：https://blog.csdn.net/victoriaw/article/details/70161180

數值方法——薄板樣條插值（Thin-Plate Spline）

插值

假設已知函式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∈R2×1a∈R2×1，向量w∈RN×1w∈RN×1，函式向量

s(x)=(σ(x−x1),σ(x−x1),⋯,σ(x−xN))Ts(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∬R2(∂2Φj∂x2)2+2(∂2Φj∂x∂y)2+(∂2Φj∂y2)2dxdyJ(Φ)=∑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