Mathematical theory of splines in tension

The traditional minimum-curvature criterion implies seeking a two-dimensional surface f(x,y) in region D, which corresponds to the minimum of the Laplacian power:  

$$\iint\limits_{D} \left\vert\nabla^2 f(x,y)\right\vert^2\,dx\,dy\;,$$ (91)

where $\nabla^2$ denotes the Laplacian operator: $\nabla^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$.

Alternatively, we can seek f(x,y) as the solution of the biharmonic differential equation  

$$(\nabla^2)^2 f(x,y) = 0\;.$$ (92)

Equation () corresponds to the normal system of equations in the least-square optimization problem (), the Laplacian operator being $\bold{D}$, and the surface f(x,y) corresponding to the unknown model $\bold{m}$. Fung (1965) and Briggs (1974) derive equation () directly from () with the help of the variational calculus and Gauss's theorem.

Formula () approximates the strain energy of a thin elastic plate Timoshenko and Woinowsky-Krieger (1968). Taking tension into account modifies both the energy formula () and the corresponding equation (). Smith and Wessel (1990) suggest the following form of the modified equation:  

$$\left[(1-\lambda) (\nabla^2)^2 - \lambda (\nabla^2)\right] f(x,y) = 0\;,$$ (93)

where the tension parameter $\lambda$ ranges from 0 to 1. The corresponding energy functional is  

$$\iint\limits_{D} \left[(1-\lambda)\,\left\vert\nabla^2 f(x,... ...\lambda\,\left\vert\nabla f(x,y)\right\vert^2\right]\,dx\,dy\;.$$ (94)

Zero tension leads to the biharmonic equation () and corresponds to the minimum curvature construction. The case of $\lambda=1$ corresponds to infinite tension. Although infinite tension is physically impossible, the resulting Laplace equation does have the physical interpretation of a steady-state temperature distribution. An important property of harmonic functions (solutions of the Laplace equation) is that they cannot have local minima and maxima in the free regions. With respect to interpolation, this means that, in the case of $\lambda=1$, the interpolation surface will be constrained to have its local extrema only at the input data locations.

Norman Sleep (2000, personal communication) points out that if the tension term $\lambda \nabla^2$ is written in the form $\nabla \cdot (\lambda \nabla)$, we can follow an analogy with heat flow and electrostatics and generalize the tension parameter $\lambda$ to a local function depending on x and y. In a more general form, $\lambda$ could be a tensor allowing for an anisotropic smoothing in some predefined directions similarly to Clapp's steering-filter method Clapp et al. (1997).

To interpolate an irregular set of data values, f_k at points (x_k,y_k), we need to solve equation () under the constraint  

$$f(x_k,y_k) = f_k\;,$$ (95)

which translates to equation () in the linear operator notation. Using the results of Chapter , we can accelerate the solution by recursive filter preconditioning. If $\bold{A}$ is the discrete filter representation of the differential operator in equation (), and we can find a minimum-phase filter $\bold{D}$ whose autocorrelation is equal to $\bold{A}$, then an appropriate preconditioning operator is a recursive inverse filtering with the filter $\bold{D}$. Formulating the problem in helical coordinates Claerbout (1998a,b) enables both the spectral factorization of $\bold{A}$ and the inverse filtering with $\bold{D}$.