CurveFitting[ArrayInterpolation]  ndimensional data interpolation (table lookup)

Calling Sequence


ArrayInterpolation(xdata, ydata, xvalues, options)
ArrayInterpolation(xydata, xvalues, options)


Parameters


xdata



a list, Array, Vector, or Matrix containing the independent coordinate(s) of each of the data points, given in one of several possible forms

ydata



a list, Array, or Vector containing the dependent coordinate of each of the data points

xydata



alternate input; a list, Array, or Matrix containing both the dependent and independent coordinates of each of the data points

xvalues



a numeric value, list, Vector, or Array containing the independent coordinate(s) of one or more points whose dependent coordinate will be approximated using interpolation

options



(optional) equation(s) of the form keyword = value, where keyword is one of method, degree, endpoints, knots, uniform, verify, extrapolate, or container.





Options


–

method = spline: Perform spline interpolation. By default, natural cubic spline interpolation is used. In the 1dimensional case, is a piecewisecubic function passing through each data point in the input. In this case, = if lies in the interval , where each is a cubic polynomial such that and for all data points in the input (where ranges from to ). The coefficients of the functions are selected such that is twice continuously differentiable (), that is, and . In addition, the "natural" condition of the spline specifies that and . In the multidimensional case, is the tensor product of such spline functions, one for each dimension. Using method=spline will produce a smoother interpolant than method=cubic ( instead of ), but is more expensive to set up and more prone to numerical instability because each segment of the spline is determined globally by the positions of all other points in the dataset.


method=linear is used by default.

•

If the options degree=d and endpoints=e are given, where d is a positive integer and e is one of natural, notaknot, or periodic, then spline interpolation will be performed using the provided degree and endpoint conditions. See Spline Continuity and End Conditions for details. These options only affect the result if method=spline is used. In the multidimensional case, the same degree and endpoint conditions are used for the splines generated in each dimension. The defaults are degree=3 and endpoints=natural, in which case natural cubic spline interpolation will be performed.

•

If splines of an even degree are being used, the option knots=data forces the use of a spline function where the spline knots are positioned on the nodes. See Spline Continuity and End Conditions for details. The default method, knots=default, defines the spline knots at the midpoints of the nodes when even degree splines are used. This option has no effect when other methods are used.

•

If the option verify=false is given, then ArrayInterpolation skips the various checks it performs to ensure correctly formatted input. This can decrease the time required to solve large problems, but will prevent the function from detecting any errors in the input. If the input is improperly sorted, contains Arrays indexed from values other than 1, contains nonrectangular or C order rtables, or is otherwise formatted incorrectly, ArrayInterpolation may return incorrect results or throw an unexpected error.

•

If the option extrapolate=e is given, where e is of type extended_numeric or truefalse, then one of the following possible extrapolation methods will be used to compute if lies outside the bounding box specified by the input:

–

extrapolate = true: Perform extrapolation using the closest valid branch of the interpolating function. In the case of method=lowest and method=highest, this is not be defined for some points, in which case undefined will be returned.

–

extrapolate = false: Do not extrapolate. An error will be thrown if any point in xvalues lies outside the bounding box specified by the input.

–

extrapolate = e, where e is of type extended_numeric : Define to be e if lies outside the bounding box specified by the input. e is commonly zero or undefined.


extrapolate=true is used by default.

•

If the option container=c is given, where c is an appropriatelysized rtable, then the computation is performed inplace and the result is returned in c. c must be of the correct size and datatype to match the output of the routine. With this option, no additional memory is allocated to store the result; this is a programmerlevel feature that can be used to reduce memory usage and decrease the time spent by Maple's garbage collector. The default is container=false, in which case Maple creates and returns a new rtable containing the result.



Description


•

The ArrayInterpolation function can interpolate numeric data in dimensions, where is any positive integer.

•

The list of independent coordinates of the data points, given in xdata, can be input in a number of different ways. xdata can be:


The preferred methods minimize memory usage and execution time by avoiding unnecessary storage and verification of redundant data. In all cases, xdata must contain real values of type numeric.

•

As an alternate form of input, a single structure xydata containing all coordinates of the data points can be entered. It can be formatted in one of the following ways:


For multidimensional data, these methods are not recommended, since space is wasted storing the full grid of independent coordinates instead of a list of all the possible coordinates in each dimension. In both cases, the coordinates must be real values of type numeric.

•

The list of values to interpolate at, given in xvalues, may be input in one of the following formats:

–

for onedimensional data, a single numeric value, or a Vector, list, or 1dimensional Array of numeric values can be input. The output will be returned in a format matching the format of the input.

–

alternatively, a list of Vectors, lists, or 1dimensional Arrays can be input. The th Vector, list, or Array in the input will be assumed to contain all of the possible th coordinates of the values to interpolate at. In this case, interpolation will be performed on an by by ... by block of points, where is the length of the th Vector or Array in the input. The output will then be returned in a Vector, Matrix, list, or Array of size by by ... by .

•

If any of the data points in xvalues lie outside the rectangular bounding box specified by the input, then extrapolation will be performed to approximate their corresponding yvalues. The method by which extrapolation is performed can be controlled by using option extrapolate; see below.

•

This routine has separate numeric methods for handling hardware and software floats. The decision about which routine to use can be controlled by setting the UseHardwareFloats environment variable. If UseHardwareFloats remains unset, then hardware floats are used if and only if Digits <= evalhf(Digits), in which case all software floats in the input will be converted to hardware floats.

•

Only computations involving numeric floating point data are supported by this routine. If the input does not contain floating point data, an error will be thrown.

•

For optimal performance, all rtables in the input should be Fortran order with rectangular storage (the default). Otherwise, a conversion will take place. All rtables in the output will be Fortran order rtables with rectangular storage.

•

This function is part of the CurveFitting package, so it can be used in the short form ArrayInterpolation(..) only after executing the command with(CurveFitting). However, it can always be accessed through the long form of the command by using CurveFitting[Interpolation](..).



Examples


An introductory example. Suppose a signal is sampled several times over a given interval of time:
>

Times := [0.00,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09, 0.10]:

>

Amplitudes := [0.6, 0.0, 0.4, 0.6, 0.3, 0.1, 0.2, 0.0, 0.1, 0.3, 0.6]:

>

pointplot(Times,Amplitudes);

ArrayInterpolation to resample the data at a higher sampling frequency:
>

NewTimes := [seq(0.001*i,i=0..100)]:

>

NewAmplitudes := ArrayInterpolation(Times,Amplitudes,NewTimes):

>

pointplot(NewTimes,NewAmplitudes);

Use a cubic spline to achieve a smoother, more realistic resampling of the data:
>

NewAmplitudes := ArrayInterpolation(Times,Amplitudes,NewTimes,method=spline):

>

pointplot(NewTimes,NewAmplitudes);

Try again, using a spline that assumes the data is sampled from a periodic waveform:
>

NewAmplitudes := ArrayInterpolation(Times,Amplitudes,NewTimes,method=spline,endpoints=periodic):

>

pointplot(NewTimes,NewAmplitudes);

A 2dimensional example: a tiny grayscale image stored in a Matrix:
>

Ranges := [[seq(1..10)],[seq(1..16)]]:

>

Img := evalf(Matrix([[255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255], [255, 255, 208, 93, 22, 0, 16, 98, 210, 255, 255, 255, 255, 255, 255, 255], [255, 196, 16, 0, 0, 0, 0, 0, 19, 200, 255, 255, 54, 54, 212, 255], [255, 61, 34, 156, 231, 255, 230, 154, 26, 58, 255, 255, 61, 0, 59, 255], [255, 5, 205, 255, 255, 255, 255, 255, 196, 7, 255, 255, 230, 91, 7, 255], [255, 20, 235, 255, 255, 255, 255, 255, 235, 36, 255, 255, 255, 233, 32, 255], [255, 125, 129, 255, 255, 255, 255, 255, 133, 186, 255, 255, 217, 89, 147, 255], [255, 255, 125, 50, 91, 99, 93, 40, 22, 121, 93, 58, 10, 140, 255, 255], [255, 255, 255, 216, 131, 75, 14, 5, 7, 12, 59, 121, 214, 255, 255, 255], [255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255]])):

>

listdensityplot(Img,style=PATCHNOGRID,axes=none,range=0..255);

Upsample it to a larger image using bilinear interpolation:
>

NewRanges := [[seq(3..52)],[seq(3..82)]]/5.0:

>

NewImg := ArrayInterpolation(Ranges,Img,NewRanges,extrapolate=255):

>

listdensityplot(NewImg,style=PATCHNOGRID,axes=none,range=0..255);

Try again, using bicubic interpolation instead for a smoother fit:
>

NewImg := ArrayInterpolation(Ranges,Img,NewRanges,extrapolate=255,method=cubic):

>

listdensityplot(NewImg,style=PATCHNOGRID,axes=none,range=0..255);

A nonuniform multidimensional example. Create some 3D mesh structures to pass through a given set of points defined by a mathematical function:
>

f := (i,j) > (3sin(i))^2(3j)^2:

>

plot3d(f,0..10,0..10,axes=normal);

Define a nonuniform grid of points, and sample f over them:
>

v := Array([0,1.5,3.5,5,6,8,9]):

>

w := Array([0,3,5,6,6.5,7.5]):

>

y := Matrix(7,6,(a,b)>evalf(f(v[a],w[b]))):

Plot the data so far:
>

pointplot3d([seq(seq([v[i],w[j],y[i,j]],j=1..6),i=1..7)],axes=normal,symbol=sphere);

Create a finer mesh to interpolate over:
>

a1 := Matrix(50,50,(i,j)>i/5):

>

a2 := Matrix(50,50,(i,j)>j/5):

>

A := ArrayTools[Concatenate](3,a1,a2):

Linear interpolation produces a quick approximation to f:
>

B := ArrayInterpolation([v,w],y,A,method=linear):

>

matrixplot(Matrix(B),axes=normal);

Nearestneighbor interpolation can also be used for quick lookup purposes:
>

B := ArrayInterpolation([v,w],y,A,method=nearest):

>

matrixplot(Matrix(B),axes=normal);

Spline interpolation produces a smoother approximation to the original function f:
>

B := ArrayInterpolation([v,w],y,A,method=spline):

>

matrixplot(Matrix(B),axes=normal);

Increasing the degree of the spline approximation can increase the smoothness of the result, but results in a longer computation time, greater numerical instability, and can cause large oscillations around the edges of a data set:
>

B := ArrayInterpolation([v,w],y,A,method=spline,degree=5):

>

matrixplot(Matrix(B),axes=normal);

Finally, a large example to illustrate a few tips for increasing the speed of computations:
>

A := evalf(Vector([seq(0..50000)])):

>

B := evalf(Vector([seq(i^(2),i=0..50000)])):

>

C := evalf(Vector([seq(0..200000)]/4)):

>

time(ArrayInterpolation(A,B,C));

 (1) 
On such a large 1dimensional example, a significant portion of the execution time is spent verifying the integrity of the input data. Disabling this verification will produce a significant speedup in the execution time of the routine, but will produce incorrect results if the input is not correctly formatted or sorted:
>

time(ArrayInterpolation(A,B,C,verify=false));

 (2) 
Asserting that the data is uniform allows a faster lookup method to be used:
>

time(ArrayInterpolation(A,B,C,verify=false,uniform=true));

 (3) 
Cubic interpolation takes longer than the default, linear method:
>

time(ArrayInterpolation(A,B,C,verify=false,uniform=true,method=cubic));

 (4) 

