Thin plate energy functional

The exact thin plate energy functional (TPEF) for a function $f(x,y)$ is

$\int _{y_{0}}^{y_{1}}\int _{x_{0}}^{x_{1}}(\kappa _{1}^{2}+\kappa _{2}^{2}){\sqrt {g}}\,dx\,dy$

where $\kappa _{1}$ and $\kappa _{2}$ are the principal curvatures of the surface mapping $f$ at the point $(x,y).$ This is the surface integral of $\kappa _{1}^{2}+\kappa _{2}^{2},$ hence the ${\sqrt {g}}$ in the integrand.

Minimizing the exact thin plate energy functional would result in a system of non-linear equations. So in practice, an approximation that results in linear systems of equations is often used. The approximation is derived by assuming that the gradient of $f$ is 0. At any point where $f_{x}=f_{y}=0,$ the first fundamental form $g_{ij}$ of the surface mapping $f$ is the identity matrix and the second fundamental form $b_{ij}$ is

${\begin{pmatrix}f_{xx}&f_{xy}\\f_{xy}&f_{yy}\end{pmatrix}}$.

We can use the formula for mean curvature $H=b_{ij}g^{ij}/2$ to determine that $H=(f_{xx}+f_{yy})/2$ and the formula for Gaussian curvature $K=b/g$ (where $b$ and $g$ are the determinants of the second and first fundamental forms, respectively) to determine that $K=f_{xx}f_{yy}-(f_{xy})^{2}.$ Since $H=(k_{1}+k_{2})/2$ and $K=k_{1}k_{2},$ the integrand of the exact TPEF equals $4H^{2}-2K.$ The expressions we just computed for the mean curvature and Gaussian curvature as functions of partial derivatives of $f$ show that the integrand of the exact TPEF is

$4H^{2}-2K=(f_{xx}+f_{yy})^{2}-2(f_{xx}f_{yy}-f_{xy}^{2})=f_{xx}^{2}+2f_{xy}^{2}+f_{yy}^{2}.$

So the approximate thin plate energy functional is

$J[f]=\int _{y_{0}}^{y_{1}}\int _{x_{0}}^{x_{1}}f_{xx}^{2}+2f_{xy}^{2}+f_{yy}^{2}\,dx\,dy.$

Rotational invariance

Rotating (x,y) by theta about z-axis to (X,Y)

Original surface with point (x,y)

Rotated surface with rotated point (X,Y)

The TPEF is rotationally invariant. This means that if all the points of the surface $z(x,y)$ are rotated by an angle $\theta$ about the $z$-axis, the TPEF at each point $(x,y)$ of the surface equals the TPEF of the rotated surface at the rotated $(x,y).$ The formula for a rotation by an angle $\theta$ about the $z$-axis is

${\binom {X}{Y}}={\begin{pmatrix}\cos \theta &-\sin \theta \\\sin \theta &\cos \theta \end{pmatrix}}{\binom {x}{y}}=R{\binom {x}{y}}.$ 1

The fact that the $z$ value of the surface at $(x,y)$ equals the $z$ value of the rotated surface at the rotated $(x,y)$ is expressed mathematically by the equation

$Z(X,Y)=z(x,y)=(z\circ xy)(X,Y)$

where $xy$ is the inverse rotation, that is, $xy(X,Y)=R^{-1}(X,Y)^{\text{T}}=R^{\text{T}}(X,Y)^{\text{T}}.$ So $Z=z\circ xy$ and the chain rule implies

$Z_{i}=z_{j}R_{ij}.$ 2

In equation (2), $Z_{0}$ means $Z_{X},$ $Z_{1}$ means $Z_{Y},$ $z_{0}$ means $z_{x},$ and $z_{1}$ means $z_{y}.$ Equation (2) and all subsequent equations in this section use non-tensor summation convention, that is, sums are taken over repeated indices in a term even if both indices are subscripts. The chain rule is also needed to differentiate equation (2) since $z_{j}$ is actually the composition $z_{j}\circ xy:$

$Z_{ik}=z_{jl}R_{kl}R_{ij}$.

Swapping the index names $j$ and $k$ yields

$Z_{ij}=z_{kl}R_{jl}R_{ik}.$ 3

Expanding the sum for each pair $i,j$ yields

${\begin{array}{lcl}Z_{XX}&=&R_{00}^{2}z_{xx}+2R_{00}R_{01}z_{xy}+R_{01}^{2}z_{yy},\\Z_{XY}&=&R_{00}R_{10}z_{xx}+(R_{00}R_{11}+R_{01}R_{10})z_{xy}+R_{01}R_{11}z_{yy},\\Z_{YY}&=&R_{10}^{2}z_{xx}+2R_{10}R_{11}z_{xy}+R_{11}^{2}z_{yy}.\end{array}}$

Computing the TPEF for the rotated surface yields

${\begin{aligned}Z_{XX}^{2}+2Z_{XY}^{2}+Z_{YY}^{2}&=(R_{11}^{2}+R_{01}^{2})z_{yy}^{2}\\&+2(R_{10}R_{11}+R_{00}R_{01})^{2}z_{xx}z_{yy}\\&+2(2R_{10}^{2}R_{11}^{2}+R_{00}^{2}R_{11}^{2}+2R_{00}R_{01}R_{10}R_{11}\\&\qquad +R_{01}^{2}R_{10}^{2}+2R_{00}^{2}R_{01}^{2})z_{xy}^{2}\\&+4(R_{10}R_{11}+R_{00}R_{01})(R_{11}^{2}+R_{01}^{2})z_{xy}z_{yy}\\&+4(R_{10}^{2}+R_{00}^{2})(R_{10}R_{11}+R_{00}R_{01})z_{xx}z_{xy}\\&+(R_{10}^{2}+R_{00}^{2})z_{xx}^{2}.\\\end{aligned}}$ 4

Inserting the coefficients of the rotation matrix $R$ from equation (1) into the right-hand side of equation (4) simplifies it to $z_{xx}^{2}+2z_{xy}^{2}+z_{yy}^{2}.$

Data fitting

The approximate thin plate energy functional can be used to fit B-spline surfaces to scattered 1D data on a 2D grid (for example, digital terrain model data). Call the grid points $(x_{i},y_{i})$ for $i=1\dots N$ (with $x_{i}\in [a,b]$ and $y_{i}\in [c,d]$) and the data values $z_{i}.$ In order to fit a uniform B-spline $f(x,y)$ to the data, the equation

$\sum _{i=1}^{N}(f(x_{i},y_{i})-z_{i})^{2}+\lambda \int _{c}^{d}\int _{a}^{b}(f_{xx}^{2}+2f_{xy}^{2}+f_{yy}^{2})\,dx\,dy$ 5

(where $\lambda$ is the "smoothing parameter") is minimized. Larger values of $\lambda$ result in a smoother surface and smaller values result in a more accurate fit to the data. The following images illustrate the results of fitting a B-spline surface to some terrain data using this method.

• 
  Original terrain data
• 
  Fitted B-spline surface with large lambda and more smoothing
• 
  Fitted B-spline surface with smaller lambda and less smoothing

The thin plate smoothing spline also minimizes equation (5), but it is much more expensive to compute than a B-spline and not as smooth (it is only $C^{1}$ at the "centers" and has unbounded second derivatives there).