薄板样条插值---TPS(Thin Plate Spline)_tps插值方法逆

薄板样条插值—TPS(Thin Plate Spline)

插值

已知一系列的观测点 $({\mathbf{x}}_1,{\mathbf{y}}_1),({\mathbf{x}}_2,{\mathbf{y}}_2),\cdots ,({\mathbf{x}}_{\mathbf{n}},{\mathbf{y}}_{\mathbf{n}})$, 这一组观测值得生成函数$Y=f(X)$未知，如何通过这一组观测值拟合出一个近似于真实生成函数的表达式。拟合出来的表达式称之为插值函数。

TPS

薄板样条的插值函数形式如下：
Φ ( x ) = c + a T x + w T s ( x ) s ( x ) = ( σ ( x − x 1 ) , σ ( x − x 2 ) , ⋯   , σ ( x − x n ) ) T σ ( x ) = ∣ ∣ x ∣ ∣ x 2 log ⁡ ∣ ∣ x ∣ ∣ x 2

$$\begin{aligned} \Phi(\bold x) &= \bold {c + a^Tx + w^Ts(x)} \\ \bold {s(x)} &= \bold{(\sigma(x - x_1), \sigma(x - x_2), \cdots, \sigma(x - x_n))^T } \\ \sigma(\bold x) & = \bold {||x||_x^2\log||x||_x^2}\\ \end{aligned}$$ Φ(x)s(x)σ(x)​=c+aTx+wTs(x)=(σ(x−x1​),σ(x−x2​),⋯,σ(x−xn​))T=∣∣x∣∣x2​log∣∣x∣∣x2​​

其中, $\mathbf{c}\in {\mathbb{R}}^{1\times 1},\mathbf{a}\in {\mathbb{R}}^{\mathbf{D}\times 1},\mathbf{w}\in {\mathbb{R}}^{\mathbf{N}\times 1}$。 从中可以看出，该函数的输出值是一个标量，也就是说，如果要针对多个维度进行插值，需要求解多个插值函数。该插值函数总共存在$N+D+1$个参数。 而每个观测点都能提供一个如下的约束，共$N$个约束条件。
$$y_k=\Phi (x_k)$$
在认为添加$D+1$个约束条件：
∑ k = 1 K w k = 0 ∑ k = 1 K w k x k 1 = 0 ⋮ ∑ k = 1 K w k x k D = 0

$$\begin{aligned} \sum_{k=1}^{K}w_k &= 0 \\ \sum_{k=1}^{K}w_k {\bf x}_k^{1} &= 0 \\ \vdots \\ \sum_{k=1}^{K}w_k {\bf x}_k^{D} &= 0 \\ \end{aligned}$$ k=1∑K​wk​k=1∑K​wk​xk1​⋮k=1∑K​wk​xkD​​=0=0=0​

令：
X = [ x 1 1 x 1 2 ⋯ x 1 D x 2 1 x 2 2 ⋯ x 2 D ⋮ ⋮ ⋱ ⋮ x N 1 x N 2 ⋯ x N D ] Y = [ y 1 y 2 ⋮ y N ] S = [ σ ( x 1 − x 1 ) σ ( x 1 − x 2 ) ⋯ σ ( x 1 − x N ) σ ( x 2 − x 1 ) σ ( x 2 − x 2 ) ⋯ σ ( x 2 − x N ) ⋮ ⋮ ⋱ ⋮ σ ( x N − x 1 ) σ ( x N − x 2 ) ⋯ σ ( x N − x N ) ]

$$\begin{aligned} X& = \begin{bmatrix} \bf{x}_1^1 & \bf{x}_1^2 &\cdots & \bf{x}_1^D\\ \bf{x}_2^1 & \bf{x}_2^2 &\cdots & \bf{x}_2^D\\ \vdots & \vdots & \ddots & \vdots\\ \bf{x}_N^1 & \bf{x}_N^2 &\cdots & \bf{x}_N^D\\ \end{bmatrix} \qquad Y = \begin{bmatrix} \bf{y}_1\\ \bf{y}_2 \\ \vdots \\ \bf{y}_N\\ \end{bmatrix}\\ S &= \begin{bmatrix} \bf{\sigma(x_1 - x_1)} & \bf{\sigma(x_1 - x_2)}&\cdots & \bf{\sigma(x_1 - x_N)}\\ \bf{\sigma(x_2 - x_1)} & \bf{\sigma(x_2 - x_2)}&\cdots & \bf{\sigma(x_2 - x_N)}\\ \vdots & \vdots & \ddots & \vdots\\ \bf{\sigma(x_N - x_1)} & \bf{\sigma(x_N - x_2)}&\cdots & \bf{\sigma(x_N - x_N)}\\ \end{bmatrix} \end{aligned}$$ XS​=⎣⎢⎢⎢⎡​x11​x21​⋮xN1​​x12​x22​⋮xN2​​⋯⋯⋱⋯​x1D​x2D​⋮xND​​⎦⎥⎥⎥⎤​Y=⎣⎢⎢⎢⎡​y1​y2​⋮yN​​⎦⎥⎥⎥⎤​=⎣⎢⎢⎢⎡​σ(x1​−x1​)σ(x2​−x1​)⋮σ(xN​−x1​)​σ(x1​−x2​)σ(x2​−x2​)⋮σ(xN​−x2​)​⋯⋯⋱⋯​σ(x1​−xN​)σ(x2​−xN​)⋮σ(xN​−xN​)​⎦⎥⎥⎥⎤​​

则约束条件构成的方程组可以改写为：
[ S 1 N X 1 N T 0 0 X T 0 0 ] [ w c a ] = Γ [ w c a ] = [ Y 0 0 ]

$$\begin{bmatrix} S & 1_N & X \\ 1_N^T&0 &0 \\ X^T& 0 &0 \\ \end{bmatrix}$$ $$\begin{bmatrix} \bf w \\ \bf c \\ \bf a \end{bmatrix}$$ = \Gamma $$\begin{bmatrix} \bf w \\ \bf c \\ \bf a \end{bmatrix}$$ = $$\begin{bmatrix} \bf Y \\ 0 \\ 0 \\ \end{bmatrix}$$ ⎣⎡​S1NT​XT​1N​00​X00​⎦⎤​⎣⎡​wca​⎦⎤​=Γ⎣⎡​wca​⎦⎤​=⎣⎡​Y00​⎦⎤​
当$\Gamma$非奇异时，这个方程组有唯一解。因此可获得参数矩阵：
$$\left[\begin{array}{c}\mathbf{w}\\ \mathbf{c}\\ \mathbf{a}\end{array}\right]={\Gamma }^{-1}\left[\begin{array}{c}\mathbf{Y}\\ 0\\ 0\end{array}\right]$$

Reference
[1] 数值方法——薄板样条插值（Thin-Plate Spline）
[2] Thin plate splines 薄板样条插值个人理解
[3] 关于Thin Plate Spline （薄板样条函数）