薄板样条插值：由稀疏关键点匹配生成稠密形变场

问题背景

在图像配准中，形变场描述的是图像中每个像素或体素应该移动到哪里。对于二维图像，可以将其表示为坐标映射函数：

$$f:\mathbb{R}^2\rightarrow \mathbb{R}^2$$

给定坐标点 $p\in \mathbb{R}^{2}$，$f(p)$ 表示它变换后的目标位置。但在实际场景中，我们通常无法直接获得完整的像素级形变场，而只能得到少量几个关键点的匹配信息，例如：源图像中的关键点 $p_i$ 与目标图像中的关键点 $q_i$ 对应。那么是否存在某种算法，可以基于已知的稀疏关键点匹配，构造整幅图像的连续形变场？薄板样条（Thin Plate Spline, TPS）就是解决这一问题的经典方法。

TPS 数学形式

---

TPS 方法可以理解为一张很薄的金属板落在若干个支点上。每个支点对应一组已知的关键点匹配关系，它限制了薄板在该位置应当到达的目标位置，并在支点之间自然弯曲，形成一个连续、平滑的表面。系统达到稳定状态时，通常会趋于弯曲能量最小的形态。

对应到图像配准中，这就意味着：TPS 在满足控制点匹配的同时，尽可能构造出整体弯曲最小、过渡最自然的非刚性形变场。假设源图像与目标图像之间已知 $n$ 对关键点匹配关系。第 $i$ 对匹配由源图像中的控制点 $p_i$ 和目标图像中的对应点 $q_i$ 组成：

$$p_i=(x_i,y_i),\quad q_i=(u_i,v_i),\quad i=1,\dots,n.$$

其中，$p_i$ 可以理解为薄板上的第 $i$ 个支点位置，$q_i$ 则规定了该支点在目标图像中应当到达的位置。因此，TPS 要构造一个连续变换函数 $f$ ，该函数既能够经过所有已知支点, 即：$f(p_i)=q_i$，又能够在支点之间保持自然、平滑的弯曲形态，使弯曲能量最小。其目标函数可以表示为：

$$\tag{1} \begin{split} \min_f \iint \left(f_{xx}^2+2f_{xy}^2 +f_{yy}^2\right)dxdy, \\ s.t.\ f(p_i)=q_i, \quad i=1,\dots,n. \end{split}$$

TPS 对应的是一个最小弯曲能量问题。对该能量泛函进行变分求解，可以得到其 Euler-Lagrange 方程，即四阶双调和方程：$\Delta^2 f = 0$。进一步考虑控制点约束后，每个控制点相当于一个点源，其响应由双调和算子的 Green 函数给出。因此，TPS 标量函数最终可以写成仿射项与若干径向基函数叠加的形式：

$$\tag{3} \begin{split} f(x,y) &= \underbrace{a_0+a_1x+a_2y}_{\text{全局仿射项}} + \sum_{i=1}^{n} \underbrace{ w_i\phi\left(\left\|(x,y)-p_i\right\|\right) }_{\text{非刚性形变项}} . \end{split}$$

其中，$p_i$ 表示第 $i$ 个控制点，$w_i$ 表示对应的径向基权重。第一项 $a_0+a_1x+a_2y$ 描述整体仿射变化，第二项描述由控制点诱导的非刚性弯曲。在二维空间中，双调和算子的 Green 函数为：

$$G(r)=\frac{1}{8\pi}r^2\log r.$$

其中，$r=\left\|(x,y)-p_i\right\|$ 表示当前位置到控制点的距离。由于前面的常系数 $\frac{1}{8\pi}$ 可以吸收到权重参数 $w_i$ 中，因此 TPS 中通常直接取径向基函数为：

$$\phi(r)=r^2\log r.$$

在二维图像配准中，目标点 $q_i=(u_i,v_i)$ 包含水平和垂直两个坐标分量，因此需要分别求解两个坐标各自的 TPS 标量函数 $f_u(x,y)$ 和 $f_v(x,y)$。

$$f(x,y)=\big(f_u(x,y),f_v(x,y)\big).$$

关于其中具体的变分过程、双调和方程的推导，以及二维 Green 函数为何对应 $r^2\log r$，将留到后续文章展开；本文先从结果形式出发，说明 TPS 如何由稀疏控制点构造连续形变场。

TPS 参数求解

---

现在我们确定了 TPS 的函数形式，剩下的问题就是求解其中的系数。为简化说明，先考虑一个标量函数：

$$f(x,y)=a_0+a_1x+a_2y+\sum_{i=1}^{k}w_i\phi(\|(x,y)-p_i\|_2),$$

其中，$p_i=(x_i,y_i)$ 为第 $i$ 个控制点，$w_i$ 为径向基函数权重，$a_0,a_1,a_2$ 为仿射项系数。将所有控制点代入 TPS 函数，可得到线性方程组。记目标标量值为：

$$[z_1, z_2, \cdots, z_k ] ,$$

记 $r_{ij}=\|p_i-p_j\|_2$，TPS 的完整线性系统可写为：

$$\begin{bmatrix} \phi(r_{11}) & \phi(r_{12}) & \cdots & \phi(r_{1k}) & 1 & x_1 & y_1\\ \phi(r_{21}) & \phi(r_{22}) & \cdots & \phi(r_{2k}) & 1 & x_2 & y_2\\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots\\ \phi(r_{k1}) & \phi(r_{k2}) & \cdots & \phi(r_{kk}) & 1 & x_k & y_k\\ 1 & 1 & \cdots & 1 & 0 & 0 & 0\\ x_1 & x_2 & \cdots & x_k & 0 & 0 & 0\\ y_1 & y_2 & \cdots & y_k & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} w_1\\ w_2\\ \vdots\\ w_k\\ a_0\\ a_1\\ a_2 \end{bmatrix} = \begin{bmatrix} z_1\\ z_2\\ \vdots\\ z_k\\ 0\\ 0\\ 0 \end{bmatrix},$$

前 $k$ 行表示 TPS 函数在所有控制点处需要满足 $f(p_i)=z_i$；最后三行是 TPS 的附加约束：

$$\sum_{i=1}^{k}w_i=0,\quad \sum_{i=1}^{k}w_ix_i=0,\quad \sum_{i=1}^{k}w_iy_i=0.$$

这些约束用于消除径向基部分中的仿射自由度，使仿射项与非刚性项分离，从而得到唯一解。若加入平滑正则项，只需将左上角的核矩阵加上 $\lambda I$，即把 $\phi(r_{ij})$ 所组成的矩阵替换为 $\Phi+\lambda I$。当 $\lambda$ 较大时，非刚性权重 $w_i$ 会被抑制，TPS 会逐渐接近纯仿射变换。对于二维图像配准，目标点通常为 $q_i=(u_i,v_i)$。此时分别令 $z_i=u_i$ 和 $z_i=v_i$，各求解一次上述线性系统，即可得到水平和垂直两个方向的 TPS 变换。