View source for Thin plate energy functional

The exact thin plate energy functional (TPEF) for a function <math>f(x,y)</math> is 
:<math>\int_{y_0}^{y_1} \int_{x_0}^{x_1} (\kappa_1^2 + \kappa_2^2) \sqrt{g} \,dx \,dy</math>

where <math>\kappa_1</math> and <math>\kappa_2</math> are the [[Principal curvature|principal curvatures]] of the surface mapping <math>f</math> at the point <math>(x,y).</math><ref name=":1">{{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}}</ref><ref>{{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}}</ref>  This is the [[Surface integral|surface integral]] of <math>\kappa_1^2 + \kappa_2^2,</math> hence the <math>\sqrt{g}</math> 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.<ref name=":1" /><ref name=":2">{{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}}</ref><ref>{{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}}</ref>  The approximation is derived by assuming that the gradient of <math>f</math> is 0.  At any point where <math>f_x = f_y =0,</math> the [[First fundamental form|first fundamental form]] <math>g_{ij}</math> of the surface mapping <math>f</math> is the identity matrix and the [[Second fundamental form|second fundamental form]] <math>b_{ij}</math> is

:<math>\begin{pmatrix} f_{xx} & f_{xy} \\ f_{xy} & f_{yy} \end{pmatrix}</math>.

We can use the formula for [[Mean curvature|mean curvature]] <math>H=b_{ij}g^{ij}/2</math><ref name=":0">{{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]}}</ref> to determine that <math>H = (f_{xx}+f_{yy})/2</math> and the formula for [[Gaussian curvature]] <math>K=b/g</math><ref name=":0" /> (where <math>b</math> and <math>g</math> are the determinants of the second and first fundamental forms, respectively) to determine that <math>K=f_{xx}f_{yy} - (f_{xy})^2.</math>  Since <math>H=(k_1+k_2)/2</math> and <math>K=k_1k_2,</math><ref name=":0" /> the integrand of the exact TPEF equals <math>4H^2 - 2K.</math>  The expressions we just computed for the mean curvature and Gaussian curvature as functions of partial derivatives of <math>f</math> show that the integrand of the exact TPEF is

:<math>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.</math>

So the approximate thin plate energy functional is

:<math>J[f] = \int_{y_0}^{y_1} \int_{x_0}^{x_1} f_{xx}^2 + 2f_{xy}^2 + f_{yy}^2 \,dx \,dy.</math>

== Rotational invariance ==
[[File:Rotated_coords1.png|thumb|Rotating (x,y) by theta about z-axis to (X,Y)]]
[[File:Rotated_coords_plus_surface1_1.png|thumb|Original surface with point (x,y)]]
[[File:Rotated_coords_plus_surface2_1.png|thumb|Rotated surface with rotated point (X,Y)]]
The TPEF is rotationally invariant.  This means that if all the points of the surface <math>z(x,y)</math> are rotated by an angle <math>\theta</math> about the <math>z</math>-axis, the TPEF at each point <math>(x,y)</math> of the surface equals the TPEF of the rotated surface at the rotated <math>(x,y).</math>  The formula for a [[Rotation matrix|rotation by an angle]] <math>\theta</math> about the <math>z</math>-axis is   

{{NumBlk|:| <math>\binom{X}{Y} = \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix}\binom{x}{y} = R\binom{x}{y}.</math> |{{EquationRef|1}}}} 

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

: <math>Z(X,Y) = z(x,y) = (z\circ xy)(X,Y)</math>  

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

{{NumBlk|:| <math>Z_i = z_j R_{ij}.</math> |{{EquationRef|2}}}} 

In equation ({{EquationNote|2}}), <math>Z_0</math> means <math>Z_X,</math> <math>Z_1</math> means <math>Z_Y,</math> <math>z_0</math> means <math>z_x,</math> and <math>z_1</math> means <math>z_y.</math>  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 <math>z_j </math> is actually the composition <math>z_j \circ xy:</math>    

: <math>Z_{ik} = z_{jl}R_{kl} R_{ij}</math>.    

Swapping the index names <math>j</math> and <math>k</math> yields    

{{NumBlk|:| <math>Z_{ij} = z_{kl}R_{jl}R_{ik}.</math> |{{EquationRef|3}}}}    

Expanding the sum for each pair <math>i,j</math> yields    

: <math>\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}</math>    

Computing the TPEF for the rotated surface yields    

{{NumBlk|:| <math> 
\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}  
</math>  |{{EquationRef|4}}}}   

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

== 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).<ref>{{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}}</ref><ref name=":2" />  Call the grid points <math>(x_i,y_i)</math> for <math>i=1\dots N</math> (with <math>x_i \in [a,b]</math> and <math>y_i \in [c,d]</math>) and the data values <math>z_i.</math>  In order to fit a uniform B-spline <math>f(x,y)</math> to the data, the equation

{{NumBlk|:| <math>\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</math> |{{EquationRef|5}}}} 

(where <math>\lambda</math> is the "smoothing parameter") is minimized.  Larger values of <math>\lambda</math> 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.

<gallery mode=packed>
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
</gallery>

The [[Polyharmonic spline#Polyharmonic smoothing splines|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 <math>C^1</math> at the "centers" and has unbounded second derivatives there).

==References==
{{reflist}}

[[Category:Splines (mathematics)]]

{{Sourceattribution|Thin plate energy functional}}