薄板样条插值

Thin Plate

薄板样条插值可以为一组对应的控制点提供光滑的插值结果，通过一组控制点，可以得到一个“平面”（不一定是二维平面），穿过这组控制点，并使得“平面”的弯曲能量（bending energy)最小。

穿过一组特征点的2D平面

图示的弯曲屏幕可以通过如下公式得到：

$f(x_i,y_i) = a_1 + a_2x_i + a_3y_i + \sum \limits_{j=1} \limits^{n}w_iU(|p_i - p_j|) \tag{1}$

前三个参数$a_1, a_2, a_3$可以看作是仿射变换, 第四个变换参数是关于使得“平面”弯曲以通过给定控制点。$U(r) = r^2\log r$是径向基函数。$|p_i - p_j|$中的$p_i$和$p_j$都是控制点坐标，（1）总共有$N$个弯曲参数$w_i$和$1+D$个仿射参数，其中$D$为控制点的维度。

我们可以化简公式（1）为：$f(x_i, y_i) = L_i(W_i|a_{1,i}a_{2,i}a_{3,i})^T $, 其中$L = [U(|p_i - p_1|), U(|p_i - p_2], \ldots , U(|p_i - pj|), U(|p_i - p_n), p_i, 1 n]$。我们接着定义$P = \begin{bmatrix} 1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ 1 & x_3 & y_3 \\ & \cdots &\\ 1 & x_n & y_n \end{bmatrix}$, $K = \begin{bmatrix} U(r_{11}) & U(r_{12}) & \cdots & U(r_{1n}) \\ U(r_{21}) & U(r_{22}) & \cdots & U(r_{2n}) \\ \vdots & \vdots & \vdots & \vdots \\ U(r_{n1}) & U(r_{n2}) & \cdots & U(r_{nn})\end{bmatrix}$。其中$(x1, y1)$表示控制点的坐标（以二维为例），$r_{1n} = |p_i - p_n|$。我们可以将$P$和$K$组合得到:$L = \begin{bmatrix} K & P \\ P^T & 0\end{bmatrix}$。这时，从公式（2）就可以推广到$N$个控制点的情形：

$V = L [W|a_1a_2a_3]^T \tag{3}$

其中$L \in \mathbb R^{(n+3) \times (n+3)}$, $W \in \mathbb R^{n \times 1}$, $K \in \mathbb R^{n \times n}$, $P \in \mathbb R^{n \times 3}$, $a_i \in \mathbb R$， $V = [v_1, v_2, \ldots, v_n, 0, 0, 0]^T$, $v_i = f(x_i, y_i)$表示在控制点$(x_i, y_i)$的“高度”。

Deformation

How to calculate the TPS parameters

由于L是一个对称矩阵，$[W|a_1a_2a_3]^T = VL^{-1}$，我们可以得到一组控制点关于某个维度的TPS参数（如上图所示，我们有7个控制点，每个控制点$f(x_i, y_i)$都有其“高度”，我们便可通过TPS拟合出一个平面，通过这些高度值）。但如果是Deformaition的话（假设我们将img1“扭曲”得到img2），我们已知的是一组控制点的对应关系，即$(x_{1,1}, y_{1,1})$与$(x_{2,1}, y_{2,1})$是对应的，我们便需要两组TPS参数来将图片“扭曲”。

$[X_2, Y_2] = L\begin{bmatrix} W_x & W_y \\ a_{1,x} & a_{1,y} \\ a_{2,x} & a_{2,y} \\ a_{3,x} & a_{3,y}\end{bmatrix}$

以嘴角为例，可以看到x方向上的“高度”为0，而y方向上的“高度”却有很大的值。

How to conduct image warping

现假设我们已经求得了一组图片之间的TPS参数$c_x \in \mathbb R^{(N+3) \times 1}, c_y \in \mathbb R^{(N+3) \times 1}$分别代表x方向和y方向。我们可以根据下面这个公式求得“扭曲”之后的点的位置（其实就是公式（1））。

$\begin{bmatrix} x’ & y’\end{bmatrix} = \begin{bmatrix} U(||(x_{i,1}, y_{i,1}) - (x_{t,m}, y_{t, n})||_2) & P\end{bmatrix}^T[c_x, c_y] \tag{4}$

其中，$(x_{i,1}, y_{i,1})， i \in [1 \ldots N] $表示上面提到的img1中的控制点，而$(x_{t,m}, y_{t, n})$表示任意一组待扭曲的图片初始像素坐标，我们这里记为img3。而$||(x_{i,1}, y_{i,1}) - (x_{t,m}, y_{t, n})||_2$则表示两个点之间的欧式距离。如果img3有M个点，则$||*||_2$得到的值其维度是$N \times M$。那么就可以得到扭曲后的坐标：$\begin{bmatrix} x’ & y’\end{bmatrix} \in \mathbb R^{M \times 2}$。代码如下：

% MATLAB
[N1,N2]=size(img3); % 获取img3的高，宽
[x,y]=meshgrid(1:N2,1:N1); % 由于初始图像未扭曲，其做为整数，可以这样获取
x=x(:);y=y(:);M=length(x);
% cx: 为x方向形变参数，其维度与img1中控制点个数有关，这里就是n_good+3,n_good就是控制点个数
fx_aff=cx(n_good+1:n_good+3)'*[ones(1,M); x'; y']; %公式（1）前三个项的实现
d2=dist2(X3b,[x y]); %X3b: (n_good * 2)  d2: (n_good * N)
fx_wrp=cx(1:n_good)'*(d2.*log(d2+eps)); %fx_wrp: (1 * N)
fx=fx_aff+fx_wrp;
fy_aff=cy(n_good+1:n_good+3)'*[ones(1,M); x'; y'];
fy_wrp=cy(1:n_good)'*(d2.*log(d2+eps));
fy=fy_aff+fy_wrp;

%%
%% 接下来就是根据扭曲后的坐标【fx, fy】进行插值
%%

Reference

Manual Registration with Thin Plates