Thin plate energy functional

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

1. Greiner, Günther (1994). "Variational Design and Fairing of Spline Surfaces" (PDF). Eurographics '94. Retrieved January 3, 2016.
2. Moreton, Henry P. (1992). "Functional Optimization for Fair Surface Design" (PDF). Computer Graphics. Retrieved January 4, 2016.
3. Eck, Matthias (1996). "Automatic reconstruction of B-splines surfaces of arbitrary topological type" (PDF). Proceedings of SIGGRAPH 96, Computer Graphics Proceedings, Annual Conference Series. Retrieved January 3, 2016.
4. Halstead, Mark (1993). "Efficient, Fair Interpolation using Catmull-Clark Surfaces" (PDF). Proceedings of the 20th annual conference on Computer graphics and interactive techniques. Retrieved January 4, 2016.
5. Kreyszig, Erwin (1991). Differential Geometry. Mineola, New York: Dover. pp. 131. ISBN 0-486-66721-9.
6. Hjelle, Oyvind (2005). "Multilevel Least Squares Approximation of Scattered Data over Binary Triangulations" (PDF). Computing and Visualization in Science. Retrieved January 14, 2016.