next up previous [pdf]

Next: Finite differences and spectral Up: Application of spectral factorization: Previous: Application of spectral factorization:

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:

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

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

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

Fung (1965) and Briggs (1974) derive equation (5) directly from (4) with the help of the variational calculus and Gauss's theorem.

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

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

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

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

Zero tension leads to the biharmonic equation (5) 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 the steering-filter method (Clapp et al., 1998).

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

$\displaystyle f(x_k,y_k) = f_k\;.$ (8)

We can accelerate the solution by recursive filter preconditioning. If $ \mathbf{A}$ is the discrete filter representation of the differential operator in equation (6) and we can find a minimum-phase filter $ \mathbf{D}$ whose autocorrelation is equal to $ \mathbf{A}$ , then an appropriate preconditioning operator is a recursive inverse filtering with the filter $ \mathbf{D}$ . The preconditioned formulation of the interpolation problem takes the form of the least-squares system (Claerbout, 2002)

$\displaystyle \mathbf{K}\, \mathbf{D}^{-1} \mathbf{p} \approx \mathbf{f}_k\;,$ (9)

where $ \mathbf{f}_k$ represents the vector of known data, $ \mathbf{K}$ is the operator of selecting the known data locations, and $ \mathbf{p}$ is the preconditioned variable: $ \mathbf{p} = \mathbf{D\, f}$ . After obtaining an iterative solution of system (9), we reconstruct the model $ \mathbf{f}$ by inverse recursive filtering: $ \mathbf{f} = \mathbf{D}^{-1}\,\mathbf{p}$ . Formulating the problem in helical coordinates (Mersereau and Dudgeon, 1974; Claerbout, 1998) enables both the spectral factorization of $ \mathbf{A}$ and the inverse filtering with $ \mathbf{D}$ .


next up previous [pdf]

Next: Finite differences and spectral Up: Application of spectral factorization: Previous: Application of spectral factorization:

2014-02-15