薄板样条插值

Thin Plate

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

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

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

$f (x_{i}, y_{i}) = a_{1} + a_{2} x_{i} + a_{3} y_{i} + \sum_{j = 1}^{n} w_{i} U (| p_{i} - p_{j} |)$

前三个参数 $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}], \dots, U (| p_{i} - p j |), U (| p_{i} - p_{n}), p_{i}, 1 n]$ 。我们接着定义 $P = [\begin{matrix} 1 & x_{1} & y_{1} \\ 1 & x_{2} & y_{2} \\ 1 & x_{3} & y_{3} \\ \dots \\ 1 & x_{n} & y_{n} \end{matrix}]$ , $K = [\begin{matrix} U (r_{11}) & U (r_{12}) & \dots & U (r_{1 n}) \\ U (r_{21}) & U (r_{22}) & \dots & U (r_{2 n}) \\ ⋮ & ⋮ & ⋮ & ⋮ \\ U (r_{n 1}) & U (r_{n 2}) & \dots & U (r_{n n}) \end{matrix}]$ 。其中 $(x 1, y 1)$ 表示控制点的坐标（以二维为例）， $r_{1 n} = | p_{i} - p_{n} |$ 。我们可以将 $P$ 和 $K$ 组合得到: $L = [\begin{matrix} K & P \\ P^{T} & 0 \end{matrix}]$ 。这时，从公式（2）就可以推广到 $N$ 个控制点的情形：

$V = L [W | a_{1} a_{2} a_{3}]^{T}$

其中 $L \in R^{(n + 3) \times (n + 3)}$ , $W \in R^{n \times 1}$ , $K \in R^{n \times n}$ , $P \in R^{n \times 3}$ , $a_{i} \in R$ ， $V = [v_{1}, v_{2}, \dots, 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_{1} a_{2} a_{3}]^{T} = V L^{- 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{matrix} W_{x} & W_{y} \\ a_{1, x} & a_{1, y} \\ a_{2, x} & a_{2, y} \\ a_{3, x} & a_{3, y} \end{matrix}]$

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

How to conduct image warping

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

$[\begin{matrix} x^{'} & y^{'} \end{matrix}] = {[\begin{matrix} U (| | (x_{i, 1}, y_{i, 1}) - (x_{t, m}, y_{t, n}) | |_{2}) & P \end{matrix}]}^{T} [c_{x}, c_{y}]$

其中， $(x_{i, 1}, y_{i, 1}) ， i \in [1 \dots 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{matrix} x^{'} & y^{'} \end{matrix}] \in 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