前言
本文是笔者阅读[1]过程中，遇到了关于Thin Plate Spline[5]相关的知识，因而查找若干资料学习后得到的一些笔记，本文主要参考[2,3,4]，希望对大家有所帮助。 如有谬误，请联系指出，转载请联系作者并注明出处。
∇ \nabla ∇ 联系方式：
e-mail: FesianXu@gmail.com
QQ: 973926198
github: https://github.com/FesianXu
知乎专栏: 计算机视觉/计算机图形理论与应用
微信公众号：
薄板样条插值
薄板样条插值（Thin Plate Spline,TPS）是插值方法的一种，是常用的2D插值方法。假如给定两张图片中一些相互对应的控制点，如何将其中一个图片进行特定的形变，使得其控制点可以与另一张图片的控制点重合，如Fig 1.1所示。当然，提供插值的方法也特别的多，TPS是其中一种技术，其有着一个基本假设
如果用一个薄钢板（只是一个比喻）的形变来模拟这种2D形变，在确保所有控制点能够尽可能匹配的情况下，怎么样才能使得钢板的弯曲量最小。
几乎所有的生物有关的形变都是可以用TPS来近似，因此TPS也经常被用于脸部关键点形变等相关的应用[1]。
为了描述整个插值过程，按照我们刚才所说的，需要定义两个项，一个是拟合项
E
Φ
\mathcal{E}_{\Phi}
EΦ，测量将源点变形后距离目标点的大小；第二个是扭曲项
E
d
\mathcal{E}_{d}
Ed，测量曲面的扭曲大小。那么有总的损失函数：
E
=
E
Φ
+
λ
E
d
(1.1)
\mathcal{E} = \mathcal{E}_{\Phi}+\lambda \mathcal{E}_{d} \tag{1.1}
E=EΦ+λEd(1.1)
其中的
λ
\lambda
λ为权值系数，控制允许非刚体形变发生的程度，不同的
λ
\lambda
λ对于整个拟合效果的影响如Fig 1.2所示。
其中有：
E
Φ
=
∑
i
=
1
N
∣
∣
Φ
(
p
i
)
−
q
i
∣
∣
2
(1.2)
\mathcal{E}_{\Phi} = \sum_{i=1}^{N}||\Phi(p_i)-q_i||^2 \tag{1.2}
EΦ=i=1∑N∣∣Φ(pi)−qi∣∣2(1.2)
E d = ∫ ∫ R 2 ( ( ∂ 2 Φ ∂ x 2 ) 2 + 2 ( ∂ 2 Φ ∂ x ∂ y ) 2 + ( ∂ 2 Φ ∂ y 2 ) 2 ) 2 d x d y (1.3) \mathcal{E}_d = \int \int_{\mathbb{R}^2} \bigg( \bigg( \dfrac{\partial^2\Phi}{\partial \mathrm{x}^2} \bigg)^2 + 2 \bigg( \dfrac{\partial^2\Phi}{\partial \mathrm{x} \partial \mathrm{y}} \bigg)^2 + \bigg( \dfrac{\partial^2\Phi}{\partial \mathrm{y}^2} \bigg)^2 \bigg)^2 \mathrm{dx}\mathrm{dy} \tag{1.3} Ed=∫∫R2((∂x2∂2Φ)2+2(∂x∂y∂2Φ)2+(∂y2∂2Φ)2)2dxdy(1.3)
其中的
N
N
N为控制点的数量，式子(1.2)很容易理解，是源目标经过形变函数
Φ
\Phi
Φ之后和目标之间的距离；而式子(1.3)是曲面扭曲的能量函数，由文献[6]中给出，最小化式子(1.1)的结果，可以推导出形变函数的唯一闭式解结果为：
Φ
(
p
)
=
M
⋅
p
+
m
0
+
∑
i
=
1
N
ω
i
U
(
∣
∣
p
−
p
i
∣
∣
)
(1.4)
\Phi(p) = \mathbf{M} \cdot p + m_0+\sum_{i=1}^{N} \omega_i U(||p-p_i||) \tag{1.4}
Φ(p)=M⋅p+m0+i=1∑NωiU(∣∣p−pi∣∣)(1.4)
其中
p
p
p为曲面上的任意一个点，有
p
=
(
x
,
y
)
T
p = (x,y)^{\mathrm{T}}
p=(x,y)T，
p
i
p_i
pi是对应域的控制点，而
M
=
(
m
1
,
m
2
)
\mathbf{M} = (m_1,m_2)
M=(m1,m2)，而这里的
U
(
⋅
)
U(\cdot)
U(⋅)为径向基函数，表示某个曲面上的点的变形会受到所有控制点变形的影响（当然，不同控制点的影响程度不一样），有
U
(
x
)
=
r
2
log
r
(1.5)
U(x) = r^2\log{r} \tag{1.5}
U(x)=r2logr(1.5)
而
ω
i
\omega_i
ωi表示对不同径向基的加权。如Fig 1.3所示，如果我们假设每个控制点都对应一个高度，也就是
(
x
i
,
y
i
)
→
v
i
(x_i,y_i)\rightarrow v_i
(xi,yi)→vi，也就是说控制点是三维空间坐标系中的自变量，而其高度是因变量，那么我们可以再继续分析式子(1.4)中的第一项和第二项。
我们发现第一项其实是尝试用一个平面 y = M ⋅ p + m 0 y = \mathbf{M} \cdot p+m_0 y=M⋅p+m0去拟合所有的目标控制点，当然这个拟合肯定不够好，因此用第二项尝试在该平面的基础上去弯曲（当然是尽可能小的弯曲），从而达到更好的拟合效果，如Fig 1.3所示。此时有未知参数 M ∈ R 2 , m 0 ∈ R \mathbf{M} \in \mathbb{R}^2, m_0 \in \mathbb{R} M∈R2,m0∈R，和 ω i , i ∈ [ 1 , N ] \omega_i, i \in [1,N] ωi,i∈[1,N]，因此一共有 1 + 2 + N 1+2+N 1+2+N个参数，其中 D = 2 D = 2 D=2是维度， N N N是控制点数目。
我们为了求解形式一般化，用以下矩阵代表之前谈到的数值，有：
P
=
[
1
x
1
y
1
1
x
2
y
2
⋮
⋮
⋮
1
x
n
y
n
]
(1.6)
\mathbf{P} = \left[ \begin{matrix} 1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ \vdots & \vdots & \vdots \\ 1 & x_n & y_n \end{matrix} \right] \tag{1.6}
P=⎣⎢⎢⎢⎡11⋮1x1x2⋮xny1y2⋮yn⎦⎥⎥⎥⎤(1.6)
其中每一行代表一个控制点坐标，该矩阵称之为控制点矩阵。
Y
=
[
v
1
v
2
⋮
v
n
0
0
0
]
(1.7)
\mathbf{Y} = \left[ \begin{matrix} v_1 \\ v_2 \\ \vdots \\ v_n \\ 0 \\ 0 \\ 0 \\ \end{matrix} \right] \tag{1.7}
Y=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡v1v2⋮vn000⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤(1.7)
该矩阵称之为高度矩阵，后面三个0是为了形式统一填充的。
K
=
[
U
(
r
11
)
U
(
r
12
)
⋯
U
(
r
21
)
U
(
r
22
)
⋯
⋯
⋯
U
(
r
N
N
)
]
(1.8)
\mathbf{K} = \left[ \begin{matrix} U(r_{11}) & U(r_{12}) & \cdots \\ U(r_{21}) & U(r_{22}) & \cdots \\ \cdots & \cdots & U(r_{NN}) \end{matrix} \right] \tag{1.8}
K=⎣⎡U(r11)U(r21)⋯U(r12)U(r22)⋯⋯⋯U(rNN)⎦⎤(1.8)
其中
r
i
j
=
∣
∣
p
i
−
p
j
∣
∣
r_{ij} = ||p_{i}-p_{j}||
rij=∣∣pi−pj∣∣表示两个控制点之间的距离。令矩阵
L
\mathbf{L}
L为：
L
=
[
K
P
P
T
0
]
∈
R
(
N
+
3
)
×
(
N
+
3
)
(1.9)
\mathbf{L} = \left[ \begin{matrix} \mathbf{K} & \mathbf{P} \\ \mathbf{P}^{\mathrm{T}} & \mathbf{0} \end{matrix} \right] \in \mathbb{R}^{(N+3) \times (N+3)} \tag{1.9}
L=[KPTP0]∈R(N+3)×(N+3)(1.9)
那么由式子(1.4)和
Φ
(
p
i
)
=
v
i
\Phi(p_i)=v_i
Φ(pi)=vi，有：
Y
=
L
(
Ω
∣
m
0
,
m
1
,
m
2
)
T
(1.10)
\mathbf{Y} = \mathbf{L} (\Omega|m_0,m_1,m_2)^{\mathrm{T}} \tag{1.10}
Y=L(Ω∣m0,m1,m2)T(1.10)
其中
Ω
=
(
ω
1
,
⋯
,
ω
N
)
\Omega = (\omega_1,\cdots,\omega_N)
Ω=(ω1,⋯,ωN)。其中的后三行引入了一组对参数的约束（虽然我并不知道这组约束的含义，有了解的朋友请在评论区赐教，谢谢）：
∑
i
=
1
N
ω
i
=
0
∑
i
=
1
N
x
i
ω
i
=
0
∑
i
=
1
N
y
i
ω
i
=
0
(1.11)
\begin{aligned} \sum_{i=1}^N \omega_i &= 0 \\ \sum_{i=1}^N x_i\omega_i &= 0 \\ \sum_{i=1}^N y_i\omega_i &= 0 \end{aligned} \tag{1.11}
i=1∑Nωii=1∑Nxiωii=1∑Nyiωi=0=0=0(1.11)
那么从式子(1.10)我们有：
(
Ω
∣
m
0
,
m
1
,
m
2
)
T
=
L
−
1
Y
(1.12)
(\Omega|m_0,m_1,m_2)^{\mathrm{T}} = \mathbf{L}^{-1}\mathbf{Y} \tag{1.12}
(Ω∣m0,m1,m2)T=L−1Y(1.12)
当然也可以通过解线性方程组(1.10)得到参数组
(
Ω
∣
m
0
,
m
1
,
m
2
)
T
(\Omega|m_0,m_1,m_2)^{\mathrm{T}}
(Ω∣m0,m1,m2)T，一旦这个参数组计算得到，那么我们的插值函数
Φ
(
p
)
\Phi(p)
Φ(p)也就已知了，只要给定平面上任意一个点，就能通过插值函数将其插值到目标平面上。
变形(deformation)
这里介绍TPS的一个主要应用，对图片的控制点进行偏移，以达到通过控制点对图像进行特定形变的目的。如Fig 2.1所示，通过拉拽嘴角的控制点（即是蓝色点），使得周围的像素，比如 A A A点移动到了 A ′ A^{\prime} A′点，此时存在位移 ( Δ x , Δ y ) (\Delta x, \Delta y) (Δx,Δy)，此时我们需要插值这个位移。 当然，对应控制点之间的移动偏移是可以知道的，记为 Δ S = { ( Δ x 1 , Δ y 1 ) , ⋯ , ( Δ x N , Δ y N ) } \Delta \mathbf{S} = \{(\Delta x_1, \Delta y_1),\cdots,(\Delta x_N, \Delta y_N)\} ΔS={(Δx1,Δy1),⋯,(ΔxN,ΔyN)}，我们要根据已知的控制点偏移去插值图片上其他任意像素点的偏移。
不妨我们把这两个位移的分量隔离开来，不考虑两个维度之间的相关性，那么可以将第一章提到的“高度”
v
i
v_i
vi在这里理解成每一个位移的分量，那么我们有两个插值函数需要预测，即是：
Δ
x
(
p
)
=
Φ
(
p
)
Δ
x
Δ
y
(
p
)
=
Φ
(
p
)
Δ
y
(2.1)
\begin{aligned} \Delta x(p) &= \Phi(p)_{\Delta x} \\ \Delta y(p) &= \Phi(p)_{\Delta y} \end{aligned} \tag{2.1}
Δx(p)Δy(p)=Φ(p)Δx=Φ(p)Δy(2.1)
假如只是选定6个控制点，分别是图片的四个角落，右眼和右侧嘴角，如Fig 2.2所示，其中红色点表示移动之前的控制点，绿色点表示移动后的控制点，我们发现只是移动了右边眼睛和右边嘴角。那么计算出来的插值函数
Φ
\Phi
Φ为：
Φ
=
[
Φ
(
p
)
Δ
x
Φ
(
p
)
Δ
y
]
(2.2)
\Phi = \left[ \begin{matrix} \Phi(p)_{\Delta x} \\ \Phi(p)_{\Delta y} \end{matrix} \right] \tag{2.2}
Φ=[Φ(p)ΔxΦ(p)Δy](2.2)
其图像如Fig 2.3所示。我们发现在四个角落，因为不存在控制点的位移，因此
Δ
x
,
Δ
y
\Delta x, \Delta y
Δx,Δy的平面没有高度变化，而嘴角向上移动，因此对应嘴角的控制点的曲面上的
Δ
y
\Delta y
Δy有着较高的高度，而对应的
Δ
x
\Delta x
Δx则没有太大的高度变化。相反的，右眼部分则是
Δ
x
\Delta x
Δx有着较为明显的高度变化，而
Δ
y
\Delta y
Δy没有。
只要得到了这个 Δ x , Δ y \Delta x, \Delta y Δx,Δy方向的插值函数，给定任意一组控制点的变化，都可以对其他非控制点的像素位置进行插值。
Update 20201028
更新一个来自知乎朋友的问题，链接在：https://zhuanlan.zhihu.com/p/227857813
ID:qing3
问题：请问为什么公式1.4的结果是标量呀，而公式1.2中把公式1.4的结果当做向量来计算的
首先我们要知道，以二维图片的变形为例子，该TPS拟合函数 Φ ( p ) \Phi(\mathbf{p}) Φ(p)，是给定一个图片中的二维坐标 p = ( x , y ) T \mathbf{p} = (x,y)^{\mathrm{T}} p=(x,y)T，去拟合该处的灰度值（当然也可以是RGB通道，只不过如果是RGB通道的话需要三个不同的拟合函数），因此该拟合函数的输出值是标量，回到式子(1.2)，我发现该式子写得有点小问题，容易造成误解，我后面改成式子(a1)，其中的 v ( ⋅ ) v(\cdot) v(⋅)表示对给定的坐标取灰度值，因此也是一个标量函数。这样子就不会有误解了。建议可以看看Fig 1.3会有更直观的感受，其预测值可以看成是高度，因此才会用“薄板变形”去形容这个过程，是很形象的。
E Φ = ∑ i = 1 N ∣ ∣ Φ ( p i ) − v ( q i ) ∣ ∣ 2 (a1) \mathcal{E}_{\Phi} = \sum_{i=1}^{N}||\Phi(p_i)-v(q_i)||^2 \tag{a1} EΦ=i=1∑N∣∣Φ(pi)−v(qi)∣∣2(a1)
Reference
[1]. Siarohin, A., Lathuilière, S., Tulyakov, S., Ricci, E., & Sebe, N. (2019). First order motion model for image animation. In Advances in Neural Information Processing Systems (pp. 7137-7147).
[2]. http://profs.etsmtl.ca/hlombaert/thinplates/
[3]. https://www.jianshu.com/p/2cc189dfbcc5#fn3
[4]. https://www.cse.wustl.edu/~taoju/cse554/lectures/lect07_Deformation2.pdf
[5]. Bookstein, F. L. (1989). Principal warps: Thin-plate splines and the decomposition of deformations. IEEE Transactions on pattern analysis and machine intelligence, 11(6), 567-585.
[6]. 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
扫一扫
3068
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
