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

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

where ${\displaystyle \kappa _{1}}$ and ${\displaystyle \kappa _{2}}$ are the principal curvatures of the surface mapping ${\displaystyle f}$ at the point ${\displaystyle (x,y).}$^[1][2] This is the surface integral of ${\displaystyle \kappa _{1}^{2}+\kappa _{2}^{2},}$ hence the ${\displaystyle {\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.^[1][3][4] The approximation is derived by assuming that the gradient of ${\displaystyle f}$ is 0. At any point where ${\displaystyle f_{x}=f_{y}=0,}$ the first fundamental form ${\displaystyle g_{ij}}$ of the surface mapping ${\displaystyle f}$ is the identity matrix and the second fundamental form ${\displaystyle b_{ij}}$ is

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

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

${\displaystyle 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

${\displaystyle 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 ${\displaystyle z(x,y)}$ are rotated by an angle ${\displaystyle \theta }$ about the ${\displaystyle z}$-axis, the TPEF at each point ${\displaystyle (x,y)}$ of the surface equals the TPEF of the rotated surface at the rotated ${\displaystyle (x,y).}$ The formula for a rotation by an angle ${\displaystyle \theta }$ about the ${\displaystyle z}$-axis is

${\displaystyle {\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 ${\displaystyle z}$ value of the surface at ${\displaystyle (x,y)}$ equals the ${\displaystyle z}$ value of the rotated surface at the rotated ${\displaystyle (x,y)}$ is expressed mathematically by the equation

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

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

${\displaystyle Z_{i}=z_{j}R_{ij}.}$

(2)

In equation (2), ${\displaystyle Z_{0}}$ means ${\displaystyle Z_{X},}$ ${\displaystyle Z_{1}}$ means ${\displaystyle Z_{Y},}$ ${\displaystyle z_{0}}$ means ${\displaystyle z_{x},}$ and ${\displaystyle z_{1}}$ means ${\displaystyle 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 ${\displaystyle z_{j}}$ is actually the composition ${\displaystyle z_{j}\circ xy:}$

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

Swapping the index names ${\displaystyle j}$ and ${\displaystyle k}$ yields

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

(3)

Expanding the sum for each pair ${\displaystyle i,j}$ yields

${\displaystyle {\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

${\displaystyle {\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 ${\displaystyle R}$ from equation (1) into the right-hand side of equation (4) simplifies it to ${\displaystyle 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).^[6][3] Call the grid points ${\displaystyle (x_{i},y_{i})}$ for ${\displaystyle i=1\dots N}$ (with ${\displaystyle x_{i}\in [a,b]}$ and ${\displaystyle y_{i}\in [c,d]}$) and the data values ${\displaystyle z_{i}.}$ In order to fit a uniform B-spline ${\displaystyle f(x,y)}$ to the data, the equation

${\displaystyle \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 ${\displaystyle \lambda }$ is the "smoothing parameter") is minimized. Larger values of ${\displaystyle \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 ${\displaystyle C^{1}}$ at the "centers" and has unbounded second derivatives there).

References