BSINTERP Function
Computes a one- or two-dimensional spline interpolant.
Usage
result = BSINTERP(xdata, fdata)
result = BSINTERP(xdata, ydata, fdata)
Input Parameters
If a one-dimensional spline is desired, then the arguments xdata and fdata are required. If a two-dimensional, tensor-product spline is desired, then xdata, ydata, and fdata are required.
xdata—Array containing the abscissas in the x-direction of the interpolation problem.
ydata—Array containing the abscissas in the y-direction of the interpolation problem.
fdata—Array containing the ordinates of the interpolation problem. If a one-dimensional spline is being computed, then fdata (i) is the data value at xdata (i). If a two-dimensional spline is being computed, then fdata is a two-dimensional array, where fdata (i, j) is the data value at (xdata (i), ydata (i)).
Returned Value
result—A structure containing information that defines the one- or two-dimensional spline.
Keywords
Double—If present and nonzero, double precision is used.
XOrder—Specifies the order of the spline in the x-direction. Default: XOrder = 4, i.e., cubic splines
YOrder—Specifies the order of the spline in the y-direction. If a one-dimensional spline is being computed, then YOrder has no effect on the computations. Default: YOrder = 4, i.e., cubic splines
XKnots—Specifies the array of knots in the x-direction to be used when computing the definition of the spline. Default: knots are selected by function BSKNOTS using its defaults
YKnots—Specifies the array of knots in the y-direction to be used when computing the definition of the spline. Default: knots are selected by function BSKNOTS using its defaults
Discussion
BSINTERP is designed to compute either a one-dimensional spline interpolant or two-dimensional, tensor-product spline interpolant to input data. The decision of whether to compute the one- or two-dimensional spline is based on the number of arguments passed to the function. Keywords are provided to allow the user to specify the order of the spline and the knots used for the spline. When computing a one-dimensional spline, the available keywords are XOrder and XKnots. When computing a two-dimensional spline, the user can specify the order and knots in x-direction and/or y-direction using keywords XOrder, XKnots, YOrder, and YKnots.
Separate discussions on one- and two-dimensional splines follow.
One-dimensional B-splines
Given the data points x = xdata, f = fdata, and the number of elements (n) in xdata and fdata, the default action of BSINTERP computes a cubic (k = 4) spline interpolant s to the data using the default knot sequence generated by the BSKNOTS Function.
Optional argument XOrder allows the user to choose the order, k, of the spline interpolant; optional argument XKnots allows user specification of knots.
Function BSINTERP is based on the routine SPLINT by de Boor (1978, p. 204).
First, BSINTERP sorts the xdata vector and stores the result in x. The elements of the fdata vector are permuted appropriately and stored in f, yielding the equivalent data (xi, fi) for i = 0 to n – 1.
The following preliminary checks are performed on the data:
xi < xi + 1 i = 0, ..., n – 2
ti < ti + k i = 0, ..., n – 1
tt ≤ ti + 1 i = 0, ..., n + k – 2
The first test checks to see that the abscissas are distinct. The second and third inequalities verify that a valid knot sequence has been specified.
In order for the interpolation matrix to be nonsingular, tk – 1 ≤ xi ≤ tn is also checked for i = 0 to n – 1. This first inequality in the last check is necessary since the method used to generate the entries of the interpolation matrix requires that the k possibly nonzero B-splines at xi:
Bj – k + 1, ..., Bj where j satisfies tj ≤ xi < tj + 1
be well-defined (that is, j – k + 1 ≥ 0).
General conditions are not known for the exact behavior of the error in spline interpolation; however, if t and x are selected properly and the data points arise from the values of a smooth (for example, Ck) function f, i.e., fi = f(xi), then the error behaves in a predictable fashion. The maximum absolute error satisfies:
where the following is true:
For more information on this, see de Boor (1978, Chapter 13) and the references therein. This function can be used in place of function CSINTERP.
The returned value for this function is a structure. This structure contains all the information to determine the spline (stored as a linear combination of B-splines) that is computed by this function. For example, the following code sequence evaluates this spline at x and returns the value in y:
y = SPVALUE(x, spline)
Two-dimensional, Tensor-product B-splines
If arguments xdata, ydata, and fdata are all included in the call to function BSINTERP, the function computes a two-dimensional, tensor-product spline interpolant. The tensor‑product spline interpolant to data {(xi, yj, fij)}, where 0 ≤ i ≤ nx – 1 and 0 ≤ j ≤ ny – 1, has the form:
where kx and ky are the orders of the splines. These numbers are defaulted to 4 but can be set to any positive integer using keywords XOrder and YOrder. Likewise, tx and ty are the corresponding knot sequences (XKnots and YKnots). These values are defaulted to the knots returned by function BSKNOTS. The algorithm requires that the following is true:
tx (kx – 1) ≤ xi ≤ tx (nx) 0 ≤ i ≤ nx – 1
ty (ky – 1) ≤ yj ≤ ty (ny) 0 ≤ j ≤ ny – 1
Tensor-product spline interpolants in two dimensions can be computed quite efficiently by solving (repeatedly) two univariate interpolation problems. The computation is motivated by the following observations:
Setting:
note that for each fixed i from 0 to nx – 1, there are ny linear equations in the same number of unknowns as can be seen below:
The same matrix appears in the previous equation.
Thus, this matrix is factored only once, then the factorization to solve the nx right-hand sides is applied. Once this is done and hmi is computed, then the coefficients cnm are solved using the relation:
for m from 0 to ny – 1, which involves one factorization and ny solutions to the different right-hand sides. This ability of function BSINTERP is based on the SPLI2D routine by de Boor (1978, p. 347).
The returned value is a structure containing all the information to determine the spline (stored in B-spline format) that is computed by this function. For example, the following code sequence evaluates this spline at (x, y) and returns the value in z:
z = SPVALUE(x, y, spline)
Example 1
In this example, a one-dimensional B-spline interpolant to function values is computed, as shown in Figure 4-5: B-Spline Interpolant. Evaluations of the computed spline are then plotted along with the original data values. Since the default settings are being used, the interpolant is determined by the not-a-knot condition (see de Boor 1978).
; Define data values.
x = FINDGEN(11)/10
f = SIN(15 * x)
; Compute interpolant.
bs = BSINTERP(x, f)
bsval = SPVALUE(FINDGEN(100)/99, bs)
; Output results.
PLOT, FINDGEN(100)/99, bsval
OPLOT, x, f, Psym = 6
 |
Example 2
In this example, a two-dimensional, tensor-product B-spline interpolant to gridded data is computed as shown in Figure 4-6: Two-Dimensional B-Spline.
; Define the abscissas in the x-direction.
x = FINDGEN(5)/4
; Define the abscissas in the y-direction.
y = FINDGEN(5)/4
; Define the sample function values.
f = FLTARR(5, 5)
FOR i=0L, 4 DO $
f(i, *) = SIN(2 * x(i)) - COS(5 * y)
; Compute the spline interpolant.
bs = BSINTERP(x, y, f)
; Use SPVALUE to evaluate the computed spline.
bsval = SPVALUE(FINDGEN(20)/19, FINDGEN(20)/19, bs)
!P.Charsize = 1.5
!P.Multi = [0, 1, 2]
WINDOW, XSize = 400, YSize = 800
; Plot the original and computed surfaces in a tall window.
SURFACE, f, x, y
SURFACE, bsval, FINDGEN(20)/19, $
FINDGEN(20)/19
 |
Warning Errors
MATH_ILL_COND_INTERP_PROB—Interpolation matrix is ill-conditioned. Solution might not be accurate.
Fatal Errors
MATH_DUPLICATE_XDATA_VALUES—The xdata values must be distinct.
MATH_YDATA_NOT_INCREASING—The ydata values must be strictly increasing.
MATH_KNOT_MULTIPLICITY—Multiplicity of the knots cannot exceed the order of the spline.
MATH_KNOT_NOT_INCREASING—Knots must be nondecreasing.
MATH_KNOT_XDATA_INTERLACING—The ith smallest element of xdata (xi) must satisfy ti ≤ xi < ti + Order, where t is the knot sequence.
MATH_XDATA_TOO_LARGE—Array xdata must satisfy
xdatai ≤ tndata, for i = 1, ..., ndata.
xdatai ≤ tndata, for i = 1, ..., ndata.
MATH_XDATA_TOO_SMALL—Array xdata must satisfy
xdatai ≥ tOrder – 1, for i = 1, ..., ndata.
xdatai ≥ tOrder – 1, for i = 1, ..., ndata.
MATH_KNOT_DATA_INTERLACING—The ith smallest element of the data arrays xdata and ydata must satisfy ti ≤ datai + Order, where t is the knot sequence.
MATH_DATA_TOO_LARGE—Data arrays xdata and ydata must satisfy
datai ≤ tnum_data, for i = 1, ..., num_data.
datai ≤ tnum_data, for i = 1, ..., num_data.
MATH_DATA_TOO_SMALL—Data arrays xdata and ydata must satisfy
datai ≥ tOrder – 1, for i = 1, ..., num_data.
datai ≥ tOrder – 1, for i = 1, ..., num_data.
