ILNumerics: Spline Interpolation in C# / VB / .NET
Spline interpolation has become the quasi standard among all available interpolation methods. It is based on piecewise cubic polynomial functions with the useful additional property of adjacent piecewise functions exposing continous second derivatives at the shared edge point of neighboring bins. This gives spline interpolation excellent smoothness properties and is the reason why splines are often preferred over simple piecewise cubic polynomials nowadays. In addition to that very efficient algorithms exist for spline interpolation. In conjunction with ILNumerics optimized and parallelized implementation spline interpolation in C# and VB is feasible even for large and multi-dimensional data.
In this topic:
- Overview
- 1-D Spline Interpolation
- N-D Spline Interpolation
- Embedding 1-D Spline Lines into N-D Space
- Limitations
Overview
ILNumerics Interpolation Toolbox offers several spline interpolation options in C# / Visual Basic / .NET. The most flexible is one dimensional spline interpolation using spline interpolator objects, like SplineInterpolatorDouble. A functional interface Interpolation.spline() exists for convenience reasons. It wraps the interpolator object and offers the same configuration options. The functional interface is described below.
For 2, 3, or higher dimensional gridded data, Interpolation.splinen() allows to interpolate gridded data of any dimensionality. If the data to be computed are not arranged in a grid but rather as individual, scattered query points, the Interpolation.splinens() function is more efficient.
One dimensional Spline Interpolation
Example (quick start): Create various spline interpolations with Interpolation.spline().
The example above uses the spline() function with various boundary conditions. It compares them to the result obtained by kriging (dashed line) for the purpose of comparison.
The table below lists all spline interpolation options:
| Function / Object | Description | When to use? |
|---|---|---|
| Interpolation.spline() | Direct access to all one-dimensional spline interpolation options. | One time computation of new values from a given data set. |
|
SplineInterpolatorDouble SplineInterpolatorComplex SplineInterpolatorSingle SplineInterpolatorFComplex |
IDisposable Class for direct access to all configuration options. | Multiple sets of new query points need to be computed from the same set of data. |
| Interpolation.interp1() |
Compute not-a-knot splines with method: spline. |
No configuration of spline end derivatives needed. |
For documentation on the general 1-dimensional interpolation functions interp1() see here.
The spline interpolator objects are described in detail in the article: one dimensional interpolation objects.
In the following the Interpolation.spline() function is described.
1-D Interpolation with Interpolation.spline()
The Interpolation.spline() function offers an interface very similar to the general 1D interpolation function interp1().
The 1st overload by default produces not-a-knot splines and marks any value of Xn, laying outside of the domain of V as NaN. The second overload allows additional configuration of such values.
Specifying known Values
Example: most simple spline interpolation.
Known values are defined by the columns of the V parameter. Optionally, X defines the positions for the elements in V. If X is omitted or null, V is assumed to be equally spaced from 0..(V.S[0] - 1). Otherwise, X must have the same length as V.S[0] and contain monotonically increasing elements with the actual positions.
V may have multiple columns. In this case, X is used on all columns and spline interpolation is performed for each column individually - if possible in parallel.
Specifying new Query Points
Example: perform parallel spline interpolation on a set of two measurements, specifying new query points explicitly.
If the optional parameter Xn is omitted or null, new query points are auto - generated for the whole definition range of V by subdividing the bins of X into halfs. If Xn is provided, it is considered a vector of arbitrary values in no specific order. Note that these query point positions are not required to lay in the range of X! By default, such values are marked as NaN in the result.
Out-of-Range Values and Extrapolation
Commonly, somehow reliable interpolation results are expected for values within the measurement range, i.e. within the range defined by X. The outOfRangeValues parameter of the second spline() overload allows one to control the value assigned to new query points outside of the given range. If this parameter is set to a numeric scalar, its value is assigned to all out-of-range elements in the result. If the parameter is set to null explicitly spline extrapolation is performed on these values.
Example: extrapolate values using spline functions
Specifying Boundary Conditions
The spline interpolation function is computed based on the known values and their positions. The resulting function will always pass through the given points V at positions X. However, the solution is not unique! An unlimited number of such spline functions exist. Users must specify one more information in order to select the concrete spline interpolating function result: the 1st derivatives (slopes) at both ends of the spline.
Three individual options for such boundary conditions exist for Interpolation.spline(). For the lower bound the derivative is configured via the lbDeriv parameter. The ubDeriv parameter defines the derivative at the upper bound:
- Not-a-knot spline (default): besides the continous second derivatives, the third derivative is made continous also. This corresponds to the same piecewise cubic function being used at the first / last bin than for the second / fore-last one. Not-a-knot splines are created by providing double.PositiveInfinity or null to the corresponding parameter.
- Natural splines create a function without curvature (second derivative = 0) at the end. Set the corresponding parameter to double.NaN for natural splines.
- Explicit derivative: the slope (first derivative) at one end can be specified explicitly by setting the corresponding parameter to a numeric value.
Example: create 3 spline interpolations based on the same data, with various boundary conditions.
Note how the 3 spline interpolations are performed in a single call to spline(). This works, because we duplicated the measurement values in V three times and individual spline interpolations are performed on the 3 columns of V. The only distinction for the individual splines is done in the lbDeriv parameter. Here, a vector is provided with 3 elements. Each element corresponds to one set of known values / one column in V and specifies the boundary condition for the corresponding column.
Note further, how the spline functions here expose their local nature: The choice of boundary condition influences the whole function - over all bins. But this influence clearly decreases as we get further away from the boundary.
N-D Spline Interpolation
Spline interpolation in ILNumerics Interpolation Toolbox is capable to perform multi-dimensional interpolation on 2D, 3D or N-dimensional datasets. Similar to the 1D case, interpolation computes new query points based on 'known' values. The known values are always expected to exist in a way that forms a rectiliniear grid. However, the new query points can similarly form such a grid, or they may exist as individual, scattered points. Both options are described in the following paragraphs.
N-D Spline Interpolation in Grids
Spline interpolation of gridded query points from gridded data in $R^n$ is performed using the Interpolation.splinen() function.
This function is capable of handling any dimensionality $n$ of V. Hence, grid position vectors X for the values in V as well as Xn for the new query points are provided as cell arrays.
Just like for the 1-dimensional spline() function, an overload exists which allows to control the value assigned to out-of-range elements or to compute extrapolated values.
Note, that no boundary condition parameters exist for Interpolation.splinen(). All spline functions computed realize not-a-knot boundary conditions.
Example: spline interpolate / smooth 2D terrain data.
Data of higher dimensionality are applied similarly.
Interpolation.splinen is called from Interpolation.interpn for method = InterpolationMethod.spline.
N-D Spline Interpolation of Scattered Points
If the new query points for interpolation are not arranged in a grid but the known points V are, spline interpolation is performed at scattered positions with the Interpolation.splinens() function.
Example: Interpolate scattered points from a grid of known values.
Interpolation.splinens() for scattered spline interpolation exposes an interface very similar to splinen() for gridded query point. If possible, the coordinates for the new query points are provided as a matrix Xn, having the scattered points defined as columns.
Note, that spline interpolation performed on scattered points that way retains the outstanding smoothness properties from the general spline() and splinens() functions. However, this comes to the price of higher memory requirements. splinens() needs memory corresponding to $k$ times the size of points stored in V, with $k \approx (prod(size(V) * 4)$. For larger data sets and/or less strict smoothness requirements, one can resort to interpns() with piecewise cubic polynomial interpolation. Another considerable alternative is given by kriging interpolation.
Dimension Order for splinen() and splinens()
Both, splinen() and splinens() define dimensions in their natural order. I.e.: the first dimension specified in the first cell element of X and Xn corresponds to the columns of V. The second cell element defined corresponds to the rows. This scheme differs from the one used by the functions interp2,interp3, interp2s, interp3s, where the 1st and the 2nd dimensions are exchanged for the purpose of easier plotting and compatibility to CAS like Octave and Matlab®. splinen and splinens expect the dimensions in the same order as interpn and interpns.
Spline Interpolating 1D Paths in N-D
The Interpolation.splinepath function allows to create a smooth line going through arbitrary points in $R^n$.
Example: create a smooth line through points in 3D.
Limitations
As is valid for all interpolation methods, spline interpolation works very satisfactory if the new query points are computed from a region which is sufficiently supported by known, close-by values. The further one gets away from known values, the more increases the uncertainty of the computed results. Unlike kriging spline interpolation does not provide any error information.
While spline interpolation produces smooth interpolation results, the interpolating function is influenced by all supporting points. As a negative side effect, any NaN values in V render the solution at all positions invalid. Therefore, NaN values in V are not supported and one must prevent from them! The following options exist for data containing NaN values:
- Resort to other interpolation methods: kriging or piecewise cubic polynomial interpolation. Latter produces results with closest smoothness properties but works on a local function base. Hence, NaN values in V affect neighboring regions only.
- Substitute NaN values with valid numeric values, obtained by interpolation.