# Thin plate spline

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

{{ safesubst:#invoke:Unsubst||$N=Confusing |date=__DATE__ |$B= {{#invoke:Message box|ambox}} }}

This is a brief derivation for the closed form solutions for smoothing Thin Plate Spline. Details about these splines can be found in (Wahba, 1990).

Thin plate splines (TPS) were introduced to geometric design by Duchon (Duchon, 1976). The name thin plate spline refers to a physical analogy involving the bending of a thin sheet of metal. In the physical setting, the deflection is in the ${\displaystyle z}$ direction, orthogonal to the plane. In order to apply this idea to the problem of coordinate transformation, one interprets the lifting of the plate as a displacement of the ${\displaystyle x}$ or ${\displaystyle y}$ coordinates within the plane. In 2D cases, given a set of ${\displaystyle K}$ corresponding points, the TPS warp is described by ${\displaystyle 2(K+3)}$ parameters which include 6 global affine motion parameters and ${\displaystyle 2K}$ coefficients for correspondences of the control points. These parameters are computed by solving a linear system, in other words, TPS has closed-form solution.

Smoothing TPS is a regularized TPS. The model has a parameter ${\displaystyle \lambda }$ to control how non-rigid is allowed for the deformation.

## Radial basis function

{{#invoke:main|main}}

Given a set of control points ${\displaystyle \{w_{i},i=1,2,\ldots ,K\}}$, a radial basis function basically defines a spatial mapping which maps any location ${\displaystyle x}$ in space to a new location ${\displaystyle f(x)}$, represented by,

${\displaystyle f(x)=\sum _{i=1}^{K}c_{i}\varphi (\left\|x-w_{i}\right\|)}$

where ${\displaystyle \left\|\cdot \right\|}$ denotes the usual Euclidean norm and ${\displaystyle \{c_{i}\}}$ is a set of mapping coefficients. One possible choice for the kernel function ${\displaystyle \varphi }$ is the thin plate spline ${\displaystyle \varphi (r)=r^{2}\log r}$. It has a more global nature than the Gaussian kernel ${\displaystyle \varphi (r)=\exp(-r^{2}/\sigma ^{2})}$, which is another common function -- a small perturbation of one of the control points always affects the coefficients corresponding to all the other points as well. Note that a thin plate spline is generally understood as a function minimizing the integral of the squared second derivative, typically in two dimensions. This corresponds to the radial basis kernel ${\displaystyle \varphi (r)=r^{2}\log r}$. Other choices of radial basis kernel produce interpolation that would not normally be described as a thin plate spline. For example the Gaussian kernel corresponds to minimization of an infinite sum of derivative terms.

## Thin plate spline

### Smoothness measure

One of the simplest smoothness measures is the space integral of the square of the second order derivatives of the mapping function. This leads us to the thin plate spline (TPS). The TPS fits a mapping function ${\displaystyle f(x)}$ between corresponding point-sets ${\displaystyle \{y_{i}\}}$ and ${\displaystyle \{x_{i}\}}$ by minimizing the following energy function:

${\displaystyle E=\iint \left[\left({\frac {\partial ^{2}f}{\partial x^{2}}}\right)^{2}+2\left({\frac {\partial ^{2}f}{\partial xy}}\right)^{2}+\left({\frac {\partial ^{2}f}{\partial y^{2}}}\right)^{2}\right]{\textrm {d}}x\,{\textrm {d}}y}$

And for a smoothing TPS, it is

${\displaystyle E_{tps}=\sum _{i=1}^{K}\|y_{i}-f(x_{i})\|^{2}+\lambda \iint \left[\left({\frac {\partial ^{2}f}{\partial x^{2}}}\right)^{2}+2\left({\frac {\partial ^{2}f}{\partial x\partial y}}\right)^{2}+\left({\frac {\partial ^{2}f}{\partial y^{2}}}\right)^{2}\right]{\textrm {d}}x\,{\textrm {d}}y}$

Then smoothing TPS is defined as

${\displaystyle f_{tps}=\arg \min _{f}E_{tps}}$

For this variational problem, it can be shown that there exists a unique minimizer ${\displaystyle f}$ (Wahba,1990) with a fixed weight parameter ${\displaystyle \lambda }$ which is presented in the next section. The finite element discretization of this variational problem, the method of elastic maps, is used for data mining and nonlinear dimensionality reduction.

### Spline

Suppose the points are in 2D (D = 2). One can use homogeneous coordinates for the point-set where a point ${\displaystyle y_{i}}$ is represented as a vector ${\displaystyle (1,y_{ix},y_{iy})}$. The unique minimizer ${\displaystyle f}$ is parameterized by ${\displaystyle \alpha }$ which comprises two matrices ${\displaystyle d}$ and ${\displaystyle c}$ (${\displaystyle \alpha =\{d,c\}}$).

${\displaystyle f_{tps}(z,\alpha )=f_{tps}(z,d,c)=z\cdot d+\sum _{i=1}^{K}\phi (\|z-x_{i}\|)\cdot c_{i}}$

where d is a ${\displaystyle (D+1)\times (D+1)}$ matrix representing the affine transformation (hence ${\displaystyle z}$ is a ${\displaystyle 1\times (D+1)}$ vector) and c is a ${\displaystyle K\times (D+1)}$ warping coefficient matrix representing the non-affine deformation. The kernel function ${\displaystyle \phi (z)}$ is a ${\displaystyle 1\times K}$ vector for each point ${\displaystyle z}$, where each entry ${\displaystyle \phi _{i}(z)=\|z-x_{i}\|^{2}\log \|z-x_{i}\|}$ for each (${\displaystyle D}$) dimensions. Note that for TPS, the control points ${\displaystyle \{w_{i}\}}$ are chosen to be the same as the set of points to be warped ${\displaystyle \{x_{i}\}}$, so we already use ${\displaystyle \{x_{i}\}}$ in the place of the control points.

If one substitutes the solution for ${\displaystyle f}$, ${\displaystyle E_{tps}}$ becomes:

${\displaystyle E_{tps}(d,c)=\|Y-Xd-\Phi c\|^{2}+\lambda {\textrm {Tr}}(c^{T}\Phi c)}$

where ${\displaystyle Y}$ and ${\displaystyle X}$ are just concatenated versions of the point coordinates ${\displaystyle y_{i}}$ and ${\displaystyle x_{i}}$, and ${\displaystyle \Phi }$ is a ${\displaystyle (K\times K)}$ matrix formed from the ${\displaystyle \phi (\|x_{i}-x_{j}\|)}$. Each row of each newly formed matrix comes from one of the original vectors. The matrix ${\displaystyle \Phi }$ represents the TPS kernel. Loosely speaking, the TPS kernel contains the information about the point-set's internal structural relationships. When it is combined with the warping coefficients ${\displaystyle c}$, a non-rigid warping is generated.

A nice property of the TPS is that it can always be decomposed into a global affine and a local non-affine component. Consequently, the TPS smoothness term is solely dependent on the non-affine components. This is a desirable property, especially when compared to other splines, since the global pose parameters included in the affine transformation are not penalized.

### Solution

The separation of the affine and non-affine warping space is done through a QR decomposition (Wahba,1990).

${\displaystyle X=[Q_{1}|Q_{2}]\left({\begin{array}{cc}R\\0\end{array}}\right)}$

where Q1 and Q2 are ${\displaystyle K\times (D+1)}$ and ${\displaystyle K\times (K-D-1)}$ orthonormal matrices, respectively. The matrix ${\displaystyle R}$ is upper triangular. With the QR decomposition in place, we have

${\displaystyle E_{tps}(\gamma ,d)=\|Q_{2}^{T}Y-Q_{2}^{T}\Phi Q_{2}\gamma \|^{2}+\|Q_{1}^{T}Y-Rd-Q_{1}^{T}\Phi Q_{2}\gamma \|^{2}+\lambda {\textrm {trace}}(\gamma ^{T}Q_{2}^{T}\Phi Q_{2}\gamma )}$

where ${\displaystyle \gamma }$ is a ${\displaystyle (K-D-1)\times (D+1)}$ matrix. Setting ${\displaystyle c=Q_{2}\gamma }$ (which in turn implies that ${\displaystyle X^{T}c=0}$) enables us to cleanly separate the first term in last third equation into a non-affine term and an affine term (first and second terms last equation respectively).

The least-squares energy function in the last equation can be first minimized w.r.t ${\displaystyle \gamma }$ and then w.r.t. ${\displaystyle d}$. By applying Tikhonov regularization we have

${\displaystyle {\hat {c}}=Q_{2}(Q_{2}^{T}\Phi Q_{2}+\lambda I_{(k-D-1)})^{-1}Q_{2}^{T}Y}$
${\displaystyle {\hat {d}}=R^{-1}Q_{1}^{T}(Y-\Phi {\hat {c}})}$

The minimum value of the TPS energy function obtained at the optimum ${\displaystyle ({\hat {c}},{\hat {d}})}$ is

${\displaystyle E_{bending}=\lambda \,{\textrm {trace}}[Q_{2}(Q_{2}^{T}\Phi Q_{2}+\lambda I_{(k-D-1)})^{-1}Q_{2}^{T}YY^{T}]}$

## Application

TPS has been widely used as the non-rigid transformation model in image alignment and shape matching.

The popularity of TPS comes from a number of advantages:

1. The interpolation is smooth with derivatives of any order.
2. The model has no free parameters that need manual tuning.
3. It has closed-form solutions for both warping and parameter estimation.
4. There is a physical explanation for its energy function.

## References

• Haili Chui: Non-Rigid Point Matching: Algorithms, Extensions and Applications. PhD Thesis, Yale University, May 2001.
• G. Wahba, 1990, Spline models for observational data. Philadelphia: Society for Industrial and Applied Mathematics.
• J. Duchon, 1976, Splines minimizing rotation invariant semi-norms in Sobolev spaces. pp 85-100, In: Constructive Theory of Functions of Several Variables, Oberwolfach 1976, W. Schempp and K. Zeller, eds., Lecture Notes in Math., Vol. 571, Springer, Berlin, 1977