Thin Plate Spline formulas
Given a set C of p 3D control points....
c
i1
=
x
i
c
i2
=
y
i
c
i3
=
z
i
,
i
∈
[
1
p
]
≡
C
p
×
3
=
x
1
y
1
z
1
x
2
y
2
z
2
⋮
x
p
y
p
z
p
...and regularization parameter λ, solve the unknown TPS weights w and a in linear equation system...
K
P
P
T
O
⋅
w
?
a
?
=
v
o
≡
L
p
3
×
p
3
⋅
x
?
p
3
×
1
=
b
p
3
×
1
..., where K, P and O are submatrices and w, a, v and o are subvectors, given by:
K
ij
=
U
c
i1
c
i2
−
c
j1
c
j2
I
ij
⋅
2
⋅
∥
i
,
j
∈
[
1
p
]
∧
≥
0
U
r
=
{
r
2
⋅
log
r
,
r
0
0,
r
=
0
=
1
p²
∑
i
=
1
p
∑
j
=
1
p
c
i1
c
i2
−
c
j1
c
j2
P
p
×
3
=
1
c
11
c
12
1
c
21
c
22
1
⋮
⋮
1
c
p1
c
p2
,
O
3
×
3
=
0
0
0
0
0
0
0
0
0
P
T
ij
=
P
ji
∥
i
∈
[
1
p
]
∧
j
∈
[
1
3
]
v
p
×
1
=
c
13
c
23
⋮
c
p3
,
o
3
×
1
=
0
0
0
w
?
p
×
1
=
w
1
w
2
⋮
w
p
,
a
?
3
×
1
=
a
1
a
2
a
3
,
Note that L, and thus also its submatrix K, is symmetric so you can calculate elements for upper triangle only and copy them to the lower one. Also,
(mean of distances between control points' xy-projections) is a constant only present on the diagonal of K, so you can easily calculate it while filling up the uppert and lower triangles.
Then, once you know values for w and a, you can interpolate z for arbitrary points (x,y) with :
z
x
,
y
=
a
1
a
2
x
a
3
y
∑
i
=
1
p
w
i
U
c
i1
c
i2
−
x
y
Function U, used both in solving the weights and interpolating, is known as the 'thin plate spline radial basis function'. Bending energy of a TPS is given by:
I
f
=
w
T
K
w
"Thin Plate Spline formulas"
newline
"Given a set C of p 3D control points.... "
newline
left lbrace stack{c_i1 = x_i # c_i2 = y_i # c_i3 = z_i } right rbrace, i in [1dotslow p] ~equiv~ C_{p times 3} = left [ matrix{x_1 # y_1 # z_1 ## x_2 # y_2 # z_2 ## ` # dotsvert # ` ## x_p # y_p # z_p } right]
newline
"...and regularization parameter λ, solve the unknown TPS weights w and a in linear equation system... "
newline
left [ matrix{ K # P ## P^T # O } right ] cdot left [ matrix{ {vec w}^? ## {vec a}^? } right ] =left [ matrix{ vec v ## vec o } right ]
~equiv~ L_{(p+3) times (p+3)} cdot {{vec x}^?}_{(p+3) times 1}={vec b}_{(p+3) times 1}
newline
"..., where K, P and O are submatrices and w, a, v and o are subvectors, given by: "
newline
K_{ij} = U left (abs{`left[ c_{i1} ` c_{i2} right] - left[ c_{j1} ` c_{j2} right] ` } right ) + I_{ij} cdot %alpha^2 cdot %lambda ~parallel ~
{i,j} in [1 dotslow p] `and` %lambda >= 0
newline
U(r)=left lbrace stack{r^2 cdot log r, r>0 #0, r=0 } right none
newline
%alpha = 1 over p² sum from{i=1} to{p} sum from{j=1} to{p} abs{`left[ c_{i1} ` c_{i2} right] - left[ c_{j1} ` c_{j2} right] ` }
newline
P_{p times 3} = left[ matrix{ 1 # c_{11} # c_{12} ## 1 # c_{21} # c_{22} ## 1 # dotsvert # dotsvert ## 1 # c_{p1} # c_{p2}} right] ,
O_{3 times 3}=left[ matrix{ 0#0#0##0#0#0##0#0#0 } right ]
newline
{P^T}_{ij} = P_{ji} ~parallel ~ i in [1 dotslow p] `and` j in [1 dotslow 3]
newline
{vec v}_{p times 1}=left[ stack{ c_{13} # c_{23} # dotsvert # c_{p3} } right ],
{vec o}_{3 times 1}=left[ stack{ 0#0#0 } right]
{{vec w}^?}_{p times 1}=left[ stack{ w_1 # w_2 # dotsvert # w_p } right],
{{vec a}^?}_{3 times 1}=left[ stack{ a_1 # a_2 # a_3 } right],
newline
"Note that L, and thus also its submatrix K, is symmetric so you can calculate elements for upper triangle only and copy them to the lower one. Also, "%alpha" (mean of distances between control points' xy-projections) is a constant only present on the diagonal of K, so you can easily calculate it while filling up the uppert and lower triangles. "
newline
"Then, once you know values for w and a, you can interpolate z for arbitrary points (x,y) with :"
newline
{hat z}(x, y) = a_{1} + a_{2}x + a_{3}y + sum from{i=1} to{p} w_{i}
U left (abs{`left[ c_{i1} ` c_{i2} right] - left[ x ` y right] ` } right )
newline
"Function U, used both in solving the weights and interpolating, is known as the 'thin plate spline radial basis function'. Bending energy of a TPS is given by: "
newline
I_f = {vec w}^T K vec w
newline