Thin Plate Spline Interpolation

Last modified: May 24, 2025

This article is written in: 🇺🇸

Thin Plate Spline Interpolation

Thin Plate Spline (TPS) interpolation is a non‑parametric, spline‑based technique for fitting a smooth surface through scattered data in two or more spatial dimensions. In its classical 2‑D form one seeks a function $f : R^{2} \to R$ that passes through specified data points while minimising the thin‑plate bending energy—the amount a thin metal sheet would bend if it were pinned at those points. The construction extends naturally to higher dimensions and higher‑order splines.

While polynomials or other radial‑basis interpolants can achieve the same pointwise accuracy, TPS is unique in that it yields the minimum possible bending energy among all twice‑differentiable functions matching the data, so the fitted surface stays as flat (i.e. as smooth) as the constraints allow. This makes TPS a staple in image warping, geometric modelling, and shape deformation tasks.

Conceptual Illustration

Imagine you have a set of control points $(x_{i}, y_{i}, z_{i})$ in 3D space, where $(x_{i}, y_{i})$ represent spatial coordinates and $z_{i}$ is the function value at that location. Thin plate spline interpolation finds a surface $z = f (x, y)$ that exactly passes through all these points. If you imagine the surface as a thin metal sheet pinned at these points, the TPS solution is the shape the sheet would naturally take if it were free to bend but not stretch, minimizing the total bending energy:

output

The resulting surface is smooth, continuous in its derivatives, and tends to flatten out smoothly between data points.

Mathematical Formulation

Given a set of $N$ data points ${(x_{i}, y_{i}, z_{i})}_{i = 1}^{N}$ , where no two points coincide, we want to find a function:

$f (x, y) = α_{0} + α_{1} x + α_{2} y + \sum_{i = 1}^{N} w_{i} ϕ (‖ (x, y) - (x_{i}, y_{i}) ‖)$

that interpolates the given data. Here:

The $α_{0}, α_{1}, α_{2}$ terms represent a polynomial of degree 1 (a plane) that gives the global trend.
The function $ϕ (r)$ is a radial basis function chosen as:

$ϕ (r) = r^{2} \ln (r)$

which is the fundamental solution associated with the thin plate spline bending energy in 2D.

The $w_{i}$ are the coefficients for the radial basis part.

This $f (x, y)$ must satisfy the interpolation conditions:

$f (x_{i}, y_{i}) = z_{i}, i = 1, \dots, N .$

Additionally, to ensure a unique solution and remove degeneracies, $f (x, y)$ must satisfy:

$\sum_{i = 1}^{N} w_{i} = \sum_{i = 1}^{N} w_{i} x_{i} = \sum_{i = 1}^{N} w_{i} y_{i} = 0$

This leads to a linear system for the unknown parameters $α_{0}, α_{1}, α_{2}, w_{1}, \dots, w_{N}$

Derivation

I. Energy Minimization:

Thin plate splines arise from minimizing a bending energy functional:

$J [f] = \int \int {(\frac{\partial^{2} f}{\partial x^{2}})}^{2} + 2 {(\frac{\partial^{2} f}{\partial x \partial y})}^{2} + {(\frac{\partial^{2} f}{\partial y^{2}})}^{2} d x d y,$ subject to the interpolation constraints $f (x_{i}, y_{i}) = z_{i}$ .

II. Variational Problem:

Solving the Euler-Lagrange equations associated with the energy minimization under the interpolation conditions yields the form of the TPS. The solution can be shown to be a polynomial of degree at most 1 plus a weighted sum of radial basis functions $ϕ (r) = r^{2} \ln (r)$ .

III. Linear System:

Substitute $f (x, y)$ into the interpolation conditions. This produces a system of $N + 3$ linear equations (for $w_{i}, α_{0}, α_{1}, α_{2}$ ):

$[\begin{matrix} 0 & P^{⊤} P & K \end{matrix}] [\begin{matrix} α w \end{matrix}] = [\begin{matrix} 0 z \end{matrix}]$

where:

$P$ is the $N \times 3$ matrix with rows $[1, x_{i}, y_{i}]$ .
$K$ is the $N \times N$ matrix with entries $K_{i j} = ϕ (‖ (x_{i}, y_{i}) - (x_{j}, y_{j}) ‖)$ .
$z$ is the vector of observed $z_{i}$ .
$α$ is $[α_{0}, α_{1}, α_{2}]^{⊤}$ and $w = [w_{1}, \dots, w_{N}]^{⊤} .$

Solving this system yields the TPS coefficients.

Algorithm (matrix form)

I. Kernel matrix

$K \in R^{N \times N}, K_{i j} = ϕ! (‖ x_{i} - x_{j} ‖_{2}) .$

II. Polynomial matrix

$P = [\begin{matrix} 1 & x_{1} & y_{1} \\ ⋮ & ⋮ & ⋮ \\ 1 & x_{N} & y_{N} \end{matrix}] \in R^{N \times 3} .$

III. Augmented linear system

$\underset{A \in R^{(N + 3) \times (N + 3)}}{\underset{⏟}{[\begin{matrix} K & P P^{T} & 0_{3 \times 3} \end{matrix}]}} [\begin{matrix} w α \end{matrix}] = [\begin{matrix} z 0 \end{matrix}], z = (z_{1}, \dots, z_{N})^{T}$

Unknowns: $w = (w_{1}, \dots, w_{N})^{! T}$ and $α = (a_{0}, a_{1}, a_{2})^{! T}$ .
$A$ is symmetric and (for distinct points) nonsingular.

IV. Solve the linear system once. Complexity is $O (N^{3})$ with dense $K$ .

V. Evaluate at any $(x, y)$ . Cost per evaluation is $O (N)$ .

Numerical stability tip Because $ϕ (r) \sim r^{2} \ln r$ is unbounded at infinity, centre and (optionally) scale the $x_{i}$ to unit box if $N$ is large or the domain is wide.

Fully-worked example (four points)

I. Data

$(0, 0, 0), (1, 0, 1), (0, 1, 1), (1, 1, 2) (N = 4)$

II. Kernel matrix $K$

Distances $r_{i j} = ‖ x_{i} - x_{j} ‖_{2}$

$\begin{array}{ccccc} 1 & 2 & 3 & 4 \\ h l i n e 1 & 0 & 1 & 1 & \sqrt{2} \\ 2 & 1 & 0 & \sqrt{2} & 1 \\ 3 & 1 & \sqrt{2} & 0 & 1 \\ 4 & \sqrt{2} & 1 & 1 & 0 \end{array}$

Compute $ϕ (r) = r^{2} \ln r$ :

$ϕ (0) = 0$ .
$ϕ (1) = 1 \cdot \ln 1 = 0$ .
$ϕ (\sqrt{2}) = 2 \ln (\sqrt{2}) = 2 \cdot 0.69314718 = 0.69314718$ .

Thus

$K = [\begin{matrix} 0 & 0 & 0 & 0.69314718 \\ 0 & 0 & 0.69314718 & 0 \\ 0 & 0.69314718 & 0 & 0 \\ 0.69314718 & 0 & 0 & 0 \end{matrix}]$

III. Polynomial matrix $P$

$P = [\begin{matrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 1 & 1 \end{matrix}]$

IV. Augmented system

$A = [\begin{matrix} K & P P^{! T} & 0_{3 \times 3} \end{matrix}], b = [\begin{matrix} 0 1 1 2 0 0 0 \end{matrix}]$

Explicitly,

$A = [\begin{matrix} 0 & 0 & 0 & \ln 2 & 1 & 0 & 0 0 & 0 & \ln 2 & 0 & 1 & 1 & 0 0 & \ln 2 & 0 & 0 & 1 & 0 & 1 \ln 2 & 0 & 0 & 0 & 1 & 1 & 1 1 & 1 & 1 & 1 & 0 & 0 & 0 0 & 1 & 0 & 1 & 0 & 0 & 0 0 & 0 & 1 & 1 & 0 & 0 & 0 \end{matrix}]$

V. Solution

Solving

$A [\begin{matrix} w α \end{matrix}] = b$

gives

$w = [\begin{matrix} 0 0 0 0 \end{matrix}]$

$α = [\begin{matrix} a_{0} a_{1} a_{2} \end{matrix}] = [\begin{matrix} 0 1 1 \end{matrix}]$

Interpretation: the four data points lie exactly on the plane

$f (x, y) = x + y,$

so the bending-energy minimiser needs no non-linear kernel part ( $w = 0$ ).

Remarks and extensions

Topic	Notes
Singular cases	If the data are exactly coplanar and you add any numerical noise, $K$ may become rank-deficient. Add a small diagonal regulariser $λ I$ in $K$ if needed.
Complexity	Solve once in $O (N^{3})$ ; thereafter evaluations are $O (N)$ . For $N ≳ 2000$ use fast methods (partition of unity, K-D trees, or the O(N) fast-TPS of Beatson & Light).
Higher dimensions	In $d$ -D the “thin-plate” energy changes and $ϕ (r)$ becomes $r^{2} \ln r$ only for $d = 2$ . In $d = 3$ one has $ϕ (r) = r$ .
Derivative continuity	The TPS interpolant is $C^{1}$ and its second derivatives are square-integrable; ideal for smoothly warping images or terrain surfaces.

Advantages

TPS yields an infinitely differentiable surface, minimizing bending energy, and producing visually pleasing, smooth interpolants.
The method exactly passes through all given data points.
TPS works with scattered data without needing a regular grid.
It generalizes easily to higher dimensions by changing the form of $ϕ (r)$ , making it extensible.

Limitations

TPS requires solving a $(N + 3) \times (N + 3)$ linear system, which can be expensive for large $N$ .
The system matrix may become ill-conditioned with many close points, often necessitating regularization like TPS smoothing splines.
Changing or adding one point affects the entire solution, as TPS is a global method, lacking local control like piecewise methods unless combined with domain decomposition techniques.