数值方法——薄板样条插值（Thin-Plate Spline）

插值

假设已知函数 $y=f(x)$ 在 $N+1$ 个点 $x_1,x_2,\cdots,x_{N+1}$ 处的函数值 $y_1,y_2,\cdots,y_{N+1}$ ，但函数的表达式 $f(x)$ 未知，那么可以通过插值函数 $p(x)$ 来逼近未知函数 $f(x)$ ，并且 $p(x)$ 必须满足

p (x k) = y k, k = 1, 2, \dots, N + 1. (1)

$p(x_k)=y_k, k=1,2,\cdots,N+1. \tag 1$

常见的插值函数的形式有多项式函数、样条函数。

多项式函数：令 $p(x)$ 为 $N$ 次多项式函数，于是 $p(x)$ 有 $N+1$ 个参数，而由公式(1)可知这 $N+1$ 个参数满足 $N+1$ 个约束条件，所以可以求出 $p(x)$ 的表达式。
样条函数：我们知道 $N$ 阶多项式函数必然有 $N-1$ 个极值点，所以得到的插值函数摆动会比较大，这有点像机器学习中的过拟合现象，可以用样条函数来避免这个问题。这里的样条函数其实就是分段函数，表示在相邻点 $x_k$ 和 $x_{k+1}$ 之间用低阶多项式函数 $S_k(x)$ 进行插值。分段线性插值和三次样条插值都属于样条插值。

TPS

本文介绍的TPS针对的是插值问题的一种特殊情况，并且TSP插值函数的形式也比较新颖。
考虑这样一个插值问题：自变量 $\bf x$ 是2维空间中的一点，函数值 $\bf y$ 也是2维空间中的一点，并且都在笛卡尔坐标系下表示。给定 $N$ 个自变量 ${\bf x}_k$ 和对应的函数值 ${\bf y}_k$ ，求插值函数

Φ (x) = [Φ 1 (x) Φ 2 (x)] ，

$\Phi(\bf x)=\left[ \begin{matrix} \Phi_1(\bf x) \\ \Phi_2(\bf x)\end{matrix}\right]，$

使得

y k = Φ (x k) . (2)

${\bf y}_k= \Phi({\bf x}_k). \tag 2$

我们可以认为是求两个插值函数 $\Phi_1(\bf x)$ 和 $\Phi_2(\bf x)$ 。

TPS插值函数形式如下：

Φ 1 (x) = c + a T x + w T s (x) (3)

$\Phi_1(\bf x)=c + \bf a^T \bf x + {\bf w}^T {\bf s}(\bf x) \tag 3$

其中 $c$ 是标量，向量 $\bf a \in \mathbb R^{2 \times 1}$ ，向量 $\bf w\in \mathbb R^{N\times 1}$ ，函数向量

s (x) = (σ (x - x 1), σ (x - x 1), \dots, σ (x - x N)) T

$\bf s(\bf x)=(\sigma(\bf x -\bf x_1), \sigma(\bf x -\bf x_1), \cdots, \sigma(\bf x -\bf x_N))^T$

σ (x) = | | x | | 22 log | | x | | 2 .

$\sigma(\bf x)=||\bf x||^2_2\log||\bf x||_2.$

$\Phi_2(\bf x)$ 和 $\Phi_1(\bf x)$ 有一样的形式。看到这里可能会产生疑问？插值函数的形式千千万，怎么就选择公式(3)这种形式呢？我们可以把一个插值函数想象成弯曲一个薄钢板，使得它穿过给定点，这样会需要一个弯曲能量：

J (Φ) = \sum j = 1 2 \iint R 2 (\partial 2 Φ j \partial x 2) 2 + 2 (\partial 2 Φ j \partial x \partial y) 2 + (\partial 2 Φ j \partial y 2) 2 d x d y

$J(\Phi) = \sum_{j=1}^2\iint_{\mathbb R^2}\left(\frac{\partial^2 \Phi_j}{\partial x^2}\right)^2 + 2\left( \frac{\partial^2 \Phi_j}{\partial x \partial y}\right)^2 + \left( \frac{\partial^2 \Phi_j}{\partial y^2}\right)^2dxdy$

那么可以证明公式(3)是使得弯曲能量最小的插值函数。参考文献[3]中给了证明过程。

TSP插值函数 $\Phi_1$ 有 $N+3$ 个参数，而条件(2)只给出了 $N$ 个约束，我们再添加三个约束：

\sum k = 1 N w k = \sum k = 1 N x x k w k = \sum k = 1 N x y k w k = 000 (4)

$\begin{align} \sum_{k=1}^N w_k = &0 \\ \sum_{k=1}^N {\bf x}_k^x w_k =& 0 \\ \sum_{k=1}^N {\bf x}_k^y w_k = & 0 \end{align} \tag 4$

${\bf x}_k^x$ 和 ${\bf x}_k^y$ 分别表示点 $\bf x$ 的 $x$ 坐标值和 $y$ 坐标值。于是(2)和(4)可以写成

⎡ ⎣ ⎢ S 1 T N X T 1 N 00 X 00 ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ w c a ⎤ ⎦ ⎥ = ⎡ ⎣ ⎢ Y x 00 ⎤ ⎦ ⎥ (5)

$\left[\begin{matrix} S & 1_N & X \\ 1_N^T & 0 &0 \\ X^T & 0 & 0 \end{matrix}\right] \left[\begin{matrix} \bf w\\ c \\ \bf a \end{matrix}\right] = \left[\begin{matrix} Y^x \\ 0 \\ 0 \end{matrix}\right] \tag 5$
其中，

(S)ij=σ(xi−xj) $(S)_{ij}=\sigma({\bf x}_i-\bf x_j)$ ，

1N $1_N$ 表示值全为1的

N $N$ 维列向量，

X = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ x x 1 x x 2 \dots x x N x y 1 x y 2 \dots x y N ⎤ ⎦ ⎥ ⎥ ⎥ ⎥

$X = \left[\begin{matrix} {\bf x}_1^x & {\bf x}_1^y\\ {\bf x}_2^x & {\bf x}_2^y \\ \cdots & \cdots \\ {\bf x}_N^x & {\bf x}_N^y \end{matrix}\right]$

Y x = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ y x 1 y x 2 \dots y x N ⎤ ⎦ ⎥ ⎥ ⎥ ⎥

$Y^x = \left[\begin{matrix} {\bf y}_1^x \\ {\bf y}_2^x \\ \cdots \\ {\bf y}_N^x \end{matrix}\right]$

我们可以令

Γ = ⎡ ⎣ ⎢ S 1 T N X T 1 N 00 X 00 ⎤ ⎦ ⎥

$\Gamma =\left[\begin{matrix} S & 1_N & X \\ 1_N^T & 0 &0 \\ X^T & 0 & 0 \end{matrix}\right]$

那么可知当 $S$ 是非奇异矩阵时， $\Gamma$ 也是非奇异矩阵，于是参数为：

⎡ ⎣ ⎢ w c a ⎤ ⎦ ⎥ = Γ - 1 ⎡ ⎣ ⎢ Y x 00 ⎤ ⎦ ⎥ (6)

$\left[\begin{matrix} \bf w\\ c \\ \bf a \end{matrix}\right]= \Gamma^{-1}\left[\begin{matrix} Y^x \\ 0 \\ 0 \end{matrix}\right] \tag 6$

可以把 $\Phi_1$ 和 $\Phi_2$ 的参数通过一个矩阵运算计算出来：

⎡ ⎣ ⎢ w x c x a x w y c y a y ⎤ ⎦ ⎥ = Γ - 1 ⎡ ⎣ ⎢ Y x 00 Y y 00 ⎤ ⎦ ⎥ (7)

$\left[\begin{matrix} {\bf w}^x & {\bf w}^y\\ c^x &c^y\\ {\bf a}^x & {\bf a}^y \end{matrix}\right]= \Gamma^{-1}\left[\begin{matrix} Y^x & Y^y\\ 0 & 0\\ 0 &0 \end{matrix}\right] \tag 7$

我们把 $\Gamma^{-1}$ 写成下面的形式：

Γ - 1 = [Γ 11 Γ 21 Γ 12 Γ 22]

$\Gamma^{-1}=\left[ \begin{matrix} \Gamma^{11} & \Gamma^{12}\\ \Gamma^{21} & \Gamma^{22} \end{matrix}\right]$

称矩阵 $\Gamma{11}$ 为弯曲能量矩阵，其秩为 $N-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

原文链接：https://blog.csdn.net/VictoriaW/article/details/70161180