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).$ {{Cite journal|url = http://mrl.nyu.edu/~elif/thesisprop/Greiner94.pdf?origin=publication_detail|title = Variational Design and Fairing of Spline Surfaces|last = Greiner|first = Günther|date = 1994|journal = Eurographics '94|access-date = January 3, 2016}}{{Cite journal|url = http://www.cs.berkeley.edu/~sequin/PAPERS/SIGGR92_MVC_MVS.pdf|title = Functional Optimization for Fair Surface Design|last = Moreton|first = Henry P.|date = 1992|journal = Computer Graphics|access-date = January 4, 2016}} 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.{{Cite journal|url = http://www.cs.hunter.cuny.edu/~ioannis/3DP_F03/PAPERS/MODELING/HOPPE/bspline.pdf|title = Automatic reconstruction of B-splines surfaces of arbitrary topological type|last = Eck|first = Matthias|date = 1996|journal = Proceedings of SIGGRAPH 96, Computer Graphics Proceedings, Annual Conference Series|access-date = January 3, 2016}}{{Cite journal|url = http://graphics.pixar.com/people/derose/publications/FairSubdivision/paper.pdf|title = Efficient, Fair Interpolation using Catmull-Clark Surfaces|last = Halstead|first = Mark|date = 1993|journal = Proceedings of the 20th annual conference on Computer graphics and interactive techniques|access-date = January 4, 2016}} 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$ {{Cite book|title = Differential Geometry|url = https://archive.org/details/differentialgeom00krey|url-access = limited|last = Kreyszig|first = Erwin|publisher = Dover|year = 1991|isbn = 0-486-66721-9|location = Mineola, New York|pages = [https://archive.org/details/differentialgeom00krey/page/n136 131]}} 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_1k_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

File:Rotated_coords1.png

File:Rotated_coords_plus_surface1_1.png

File:Rotated_coords_plus_surface2_1.png

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

{{NumBlk2|:| $\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

{{NumBlk2|:| $Z_i = z_j R_{ij}.$ |2}}

In equation ({{EquationNote|2}}), $Z_0$ means $Z_X,$ $Z_1$ means $Z_Y,$ $z_0$ means $z_x,$ and $z_1$ means $z_y.$ Equation ({{EquationNote|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 ({{EquationNote|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

{{NumBlk2|:| $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

{{NumBlk2|:| $\begin{align}
Z_{XX}^2 + 2 Z_{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}^2R_{11}^2 + R_{00}^2R_{11}^2 + 2R_{00}R_{01}R_{10}R_{11} \\
&\qquad + R_{01}^2R_{10}^2 + 2R_{00}^2R_{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{align}$ |4}}

Inserting the coefficients of the rotation matrix $R$ from equation ({{EquationNote|1}}) into the right-hand side of equation ({{EquationNote|4}}) simplifies it to $z_{xx}^2 + 2 z_{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).{{Cite journal|url = http://home.simula.no/~oyvindhj/pub/lsbintrSpringer.pdf|title = Multilevel Least Squares Approximation of Scattered Data over Binary Triangulations|last = Hjelle|first = Oyvind|date = 2005|journal = Computing and Visualization in Science|access-date = January 14, 2016}} 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

{{NumBlk2|:| $\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.

File:Original terrain data1.png|Original terrain data

File:Fitted bspline large lambda.png|Fitted B-spline surface with large lambda and more smoothing

File:Fitted bspline smaller lambda.png|Fitted B-spline surface with smaller lambda and less smoothing

The thin plate smoothing spline also minimizes equation ({{EquationNote|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).

References

Category:Splines (mathematics)