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