Interpolation on rectilinear grids

Grid-based interpolation on rectilinear grids provides efficient evaluation of functions sampled on structured Cartesian grids in 2D, 3D, and higher dimensions. The GridSurface class and associated factory methods in the Interpolation class enable fast, accurate interpolation from gridded data.

Rectilinear Grids

A rectilinear grid (also called a product grid) is defined by independent coordinate vectors along each dimension. Grid points are located at the Cartesian product of these coordinate vectors. For example, a 2D rectilinear grid with x-coordinates [0, 1, 3] and y-coordinates [0, 2, 4, 5] has 3 × 4 = 12 grid points.

Rectilinear grids are particularly well-suited for data that naturally aligns with coordinate axes but may have non-uniform spacing. Common applications include:

Physical simulations on adaptive grids
Image or voxel data resampling
Lookup tables and tabulated functions
Geospatial elevation and gridded sensor data

Regular vs. Rectilinear Grids: A regular grid is a special case of a rectilinear grid where all axes have uniform spacing. Rectilinear grids allow non-uniform spacing along each axis, providing greater flexibility for representing data with varying resolution.

Supported Interpolation Methods

The grid interpolation API supports multiple interpolation methods, each with different trade-offs in smoothness, accuracy, and computational cost:

Grid interpolation method comparison

Method	Continuity	Dimensions	Use Cases
Nearest	Discontinuous (piecewise constant)	Any	Categorical data, extremely fast lookups
Linear	C⁰ (continuous)	Any	General-purpose, balanced performance and quality
Cubic (specialized)	C² (smooth)	2D, 3D only	High-quality visualization, natural boundary conditions
Tensor-product splines	C¹ (smooth with cubic kernel)	Any	N-D interpolation, custom kernels, flexible boundaries

Nearest-neighbor interpolation returns the value at the closest grid point. It is piecewise constant, creating a blocky appearance but offering the fastest evaluation. Use for categorical data or when minimal computational cost is essential.

Linear (multilinear) interpolation performs bilinear (2D), trilinear (3D), or multilinear (N-D) interpolation within each grid cell. The result is continuous but has discontinuous first derivatives at cell boundaries. This is the most commonly used method, providing a good balance of accuracy and speed.

Cubic spline interpolation comes in two forms:

Specialized 2D/3D cubic surfaces (created via GridCubic2DInterpolator(IReadOnlyList<Double>, IReadOnlyList<Double>, IReadOnlyList<Double>, ExtrapolationMode) and GridCubic3DInterpolator) : These precompute partial derivatives using natural cubic splines along each axis, achieving C² continuity (smooth first and second derivatives). Available only for 2D and 3D grids. Ideal for high-quality visualization and when natural boundary conditions (zero second derivative at edges) are appropriate.
Tensor-product cubic splines (created via GridCubicTensorSplineInterpolator): These use recursive 1D cubic Hermite interpolation (Catmull-Rom style) along each axis, achieving C¹ continuity. Work for any dimension including 1D, 2D, 3D, and beyond. Support custom boundary conditions (periodic, reflect, clamp, fill). Slightly less smooth than specialized 2D/3D surfaces but more flexible.

Both cubic methods can introduce overshoots and oscillations near sharp features. They are best suited for smooth underlying functions.

Specialized 2D/3D Cubic vs. Tensor-Product Splines

The library provides two families of cubic interpolation for rectilinear grids, each with distinct characteristics:

Comparison: Specialized cubic vs. tensor-product splines

Feature	Specialized 2D/3D Cubic	Tensor-Product Splines
Factory methods	GridCubic2DInterpolator(IReadOnlyList<Double>, IReadOnlyList<Double>, IReadOnlyList<Double>, ExtrapolationMode), GridCubic3DInterpolator	GridCubicTensorSplineInterpolator
Dimensions supported	2D and 3D only	Any (1D, 2D, 3D, 4D, ...)
Smoothness	C² (continuous value, 1st and 2nd derivatives)	C¹ with cubic kernel (continuous value and 1st derivative)
Underlying method	Natural cubic splines (CubicSpline) with global tridiagonal solve	Local cubic Hermite (Catmull-Rom) or custom kernels (ISplineKernel1D)
Boundary conditions	Natural only (zero 2nd derivative at boundaries)	Flexible: clamp, periodic, reflect, fill
Custom kernels	Not supported	Yes, via ISplineKernel1D
Performance	Faster evaluation (O(1) after precomputation)	Good (recursive 1D evaluation, zero allocations)

When to use specialized 2D/3D cubic surfaces:

You have 2D or 3D data
You need the highest smoothness (C²)
Natural boundary conditions are acceptable
You prioritize evaluation performance

When to use tensor-product splines:

You have N-dimensional data (any dimension)
You need custom boundary conditions (periodic, reflect, etc.)
You want to use custom spline kernels (linear, cubic, PCHIP, etc.)
C¹ continuity is sufficient

Creating Grid-Based Interpolators

The Interpolation class provides factory methods to create grid interpolators for 2D, 3D, and N-dimensional data. The methods follow a consistent naming pattern:

GridLinearInterpolator – multilinear interpolation
GridNearestInterpolator – nearest-neighbor interpolation
GridCubic2DInterpolator(IReadOnlyList<Double>, IReadOnlyList<Double>, IReadOnlyList<Double>, ExtrapolationMode), GridCubic3DInterpolator, GridCubicTensorSplineInterpolator(IReadOnlyList<IReadOnlyList<Double>>, IReadOnlyList<Double>, ExtrapolationMode) – cubic spline interpolation

Array Dimension Order

Important: When passing grid values as 2D or 3D arrays, the dimension order must match the order of the axis arguments.

For a 2D grid created with GridXxxInterpolator(xAxis, yAxis, values):

values.GetLength(0) must equal xAxis.Length
values.GetLength(1) must equal yAxis.Length
values[i, j] corresponds to the function value at f(xAxis[i], yAxis[j])

For a 3D grid, the rule extends analogously:

values[i, j, k] = f(xAxis[i], yAxis[j], zAxis[k])

This convention ensures that the first array dimension corresponds to the first axis argument, the second array dimension to the second axis, and so on. When initializing arrays, structure your loops accordingly:

double[] xAxis = { 0.0, 1.0, 2.0 };
double[] yAxis = { 0.0, 1.0 };
double[,] values = new double[xAxis.Length, yAxis.Length];

for (int i = 0; i < xAxis.Length; i++)
{
    for (int j = 0; j < yAxis.Length; j++)
    {
        values[i, j] = MyFunction(xAxis[i], yAxis[j]);
    }
}

Visual Basic

Dim xAxis() As Double = { 0.0, 1.0, 2.0 }
Dim yAxis() As Double = { 0.0, 1.0 }
Dim values(xAxis.Length - 1, yAxis.Length - 1) As Double

For i As Integer = 0 To xAxis.Length - 1
    For j As Integer = 0 To yAxis.Length - 1
        values(i, j) = MyFunction(xAxis(i), yAxis(j))
    Next
Next

Visual Basic

No code example is currently available or this language may not be supported.

let xAxis = [| 0.0; 1.0; 2.0 |]
let yAxis = [| 0.0; 1.0 |]
let values = Array2D.init xAxis.Length yAxis.Length (fun i j -> 
    myFunction xAxis.[i] yAxis.[j])

2D and 3D Examples

2D Grid Interpolation:

To create a 2D grid interpolator, provide the x-axis, y-axis, and a 2D array or matrix of function values. The values array must have dimensions matching the axis lengths as described in the Array Dimension Order section above.

// Define a 2D rectilinear grid
double[] xAxis = { 0.0, 1.0, 2.0, 3.0 };
double[] yAxis = { 0.0, 1.0, 2.0 };

// Define values on the grid
// Array dimension order: values[i, j] = f(xAxis[i], yAxis[j])
double[,] values = new double[xAxis.Length, yAxis.Length];
for (int i = 0; i < xAxis.Length; i++)
{
    for (int j = 0; j < yAxis.Length; j++)
    {
        values[i, j] = i + j * 4;  // Example function
    }
}

// Create a linear grid interpolator
var surface = Numerics.NET.Interpolation.GridLinearInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Throw);

// Evaluate at a single point
double z1 = surface.Evaluate(1.5, 1.0);
Console.WriteLine($"Linear interpolation at (1.5, 1.0): {z1}");

// Evaluate at multiple points
double z2 = surface.Evaluate(0.5, 0.5);
double z3 = surface.Evaluate(2.5, 1.5);
Console.WriteLine($"Values: ({z2:F2}, {z3:F2})");

Visual Basic

' Define a 2D rectilinear grid
Dim xAxis() As Double = { 0.0, 1.0, 2.0, 3.0 }
Dim yAxis() As Double = { 0.0, 1.0, 2.0 }

' Define values on the grid
' Array dimension order: values(i, j) = f(xAxis(i), yAxis(j))
Dim values(xAxis.Length - 1, yAxis.Length - 1) As Double
For i As Integer = 0 To xAxis.Length - 1
    For j As Integer = 0 To yAxis.Length - 1
        values(i, j) = i + j * 4  ' Example function
    Next
Next

' Create a linear grid interpolator
Dim surface = Numerics.NET.Interpolation.GridLinearInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Throw)

' Evaluate at a single point
Dim z1 = surface.Evaluate(1.5, 1.0)
Console.WriteLine($"Linear interpolation at (1.5, 1.0): {z1}")

' Evaluate at multiple points
Dim z2 = surface.Evaluate(0.5, 0.5)
Dim z3 = surface.Evaluate(2.5, 1.5)
Console.WriteLine($"Values: ({z2:F2}, {z3:F2})")

Visual Basic

No code example is currently available or this language may not be supported.

// Define a 2D rectilinear grid
let xAxis = [| 0.0; 1.0; 2.0; 3.0 |]
let yAxis = [| 0.0; 1.0; 2.0 |]

// Define values on the grid
// Array dimension order: values.[i, j] = f(xAxis.[i], yAxis.[j])
let values = Array2D.init xAxis.Length yAxis.Length (fun i j -> 
    float i + float j * 4.0)  // Example function

// Create a linear grid interpolator
let surface = Numerics.NET.Interpolation.GridLinearInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Throw)

// Evaluate at a single point
let z1 = surface.Evaluate(1.5, 1.0)
printfn $"Linear interpolation at (1.5, 1.0): {z1}"

// Evaluate at multiple points
let z2 = surface.Evaluate(0.5, 0.5)
let z3 = surface.Evaluate(2.5, 1.5)
printfn $"Values: ({z2:F2}, {z3:F2})"

For smoother interpolation, use cubic splines:

// Create a cubic spline interpolator for smoother results
var cubicSurface = Numerics.NET.Interpolation.GridCubicInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Throw);

// Evaluate at the same points - cubic will be smoother
double zc1 = cubicSurface.Evaluate(1.5, 1.0);
double zc2 = cubicSurface.Evaluate(0.5, 0.5);
Console.WriteLine($"Cubic interpolation at (1.5, 1.0): {zc1}");
Console.WriteLine($"Cubic interpolation at (0.5, 0.5): {zc2}");

Visual Basic

' Create a cubic spline interpolator for smoother results
Dim cubicSurface = Numerics.NET.Interpolation.GridCubicInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Throw)

' Evaluate at the same points - cubic will be smoother
Dim zc1 = cubicSurface.Evaluate(1.5, 1.0)
Dim zc2 = cubicSurface.Evaluate(0.5, 0.5)
Console.WriteLine($"Cubic interpolation at (1.5, 1.0): {zc1}")
Console.WriteLine($"Cubic interpolation at (0.5, 0.5): {zc2}")

Visual Basic

No code example is currently available or this language may not be supported.

// Create a cubic spline interpolator for smoother results
let cubicSurface = Numerics.NET.Interpolation.GridCubicInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Throw)

// Evaluate at the same points - cubic will be smoother
let zc1 = cubicSurface.Evaluate(1.5, 1.0)
let zc2 = cubicSurface.Evaluate(0.5, 0.5)
printfn $"Cubic interpolation at (1.5, 1.0): {zc1}"
printfn $"Cubic interpolation at (0.5, 0.5): {zc2}"

N-Dimensional Grid Interpolation:

For 3D and higher-dimensional grids, pass an array of axis coordinate vectors and a flat array of values. The values must be stored in row-major order, where the last dimension varies fastest. For a 3D grid with axes [z, y, x], the value at grid point (i, j, k) is stored at index i*ny*nx + j*nx + k.

// Create a 3D grid interpolator
Vector<double> xAxis = Vector.Create(0.0, 1.0, 2.0);
Vector<double> yAxis = Vector.Create(0.0, 1.0);
Vector<double> zAxis = Vector.Create(0.0, 1.0, 2.0);

// Axes are ordered as [z, y, x] for row-major storage
Vector<double>[] axes = new Vector<double>[] { zAxis, yAxis, xAxis };

// Flat values array in row-major order (z varies slowest, x varies fastest)
Vector<double> values = Vector.Create(
    // z=0, y=0: x values
    1.0, 2.0, 3.0,
    // z=0, y=1: x values
    4.0, 5.0, 6.0,
    // z=1, y=0: x values
    7.0, 8.0, 9.0,
    // z=1, y=1: x values
    10.0, 11.0, 12.0,
    // z=2, y=0: x values
    13.0, 14.0, 15.0,
    // z=2, y=1: x values
    16.0, 17.0, 18.0
);

// Create the ND linear interpolator
var field = Numerics.NET.Interpolation.GridLinearInterpolator(
    axes, values, ExtrapolationMode.Throw);

// Evaluate at a 3D point
Vector<double> point = Vector.Create(1.0, 0.5, 1.5);
double value = field.Evaluate(point);
Console.WriteLine($"3D interpolation at (1.0, 0.5, 1.5): {value}");

Visual Basic

' Create a 3D grid interpolator
Dim xAxis = Vector.Create(0.0, 1.0, 2.0)
Dim yAxis = Vector.Create(0.0, 1.0)
Dim zAxis = Vector.Create(0.0, 1.0, 2.0)

' Axes are ordered as [z, y, x] for row-major storage
Dim axes = New Vector(Of Double)() { zAxis, yAxis, xAxis }

' Flat values array in row-major order (z varies slowest, x varies fastest)
Dim values = Vector.Create(
    1.0, 2.0, 3.0,
    4.0, 5.0, 6.0,
    7.0, 8.0, 9.0,
    10.0, 11.0, 12.0,
    13.0, 14.0, 15.0,
    16.0, 17.0, 18.0
)

' Create the ND linear interpolator
Dim field = Numerics.NET.Interpolation.GridLinearInterpolator(
    axes, values, ExtrapolationMode.Throw)

' Evaluate at a 3D point
Dim point = Vector.Create(1.0, 0.5, 1.5)
Dim value = field.Evaluate(point)
Console.WriteLine($"3D interpolation at (1.0, 0.5, 1.5): {value}")

Visual Basic

No code example is currently available or this language may not be supported.

// Create a 3D grid interpolator
let xAxis = Vector.Create(0.0, 1.0, 2.0)
let yAxis = Vector.Create(0.0, 1.0)
let zAxis = Vector.Create(0.0, 1.0, 2.0)

// Axes are ordered as [z, y, x] for row-major storage
let axes = [| zAxis; yAxis; xAxis |]

// Flat values array in row-major order (z varies slowest, x varies fastest)
let values = Vector.Create(
    1.0, 2.0, 3.0,
    4.0, 5.0, 6.0,
    7.0, 8.0, 9.0,
    10.0, 11.0, 12.0,
    13.0, 14.0, 15.0,
    16.0, 17.0, 18.0
)

// Create the ND linear interpolator
let field = Numerics.NET.Interpolation.GridLinearInterpolator(
    axes, values, ExtrapolationMode.Throw)

// Evaluate at a 3D point
let point = Vector.Create(1.0, 0.5, 1.5)
let value = field.Evaluate(point)
printfn $"3D interpolation at (1.0, 0.5, 1.5): {value}"

Evaluating Grid Surfaces

Once a GridSurface is created, evaluation is straightforward using the Evaluate(IReadOnlyList<Double>) method.

For 2D surfaces, call surface.Evaluate(x, y). For N-dimensional surfaces, pass a coordinate array or span: surface.Evaluate(new double[] { x, y, z }).

Evaluation is efficient: linear interpolation requires a small number of arithmetic operations per query, while cubic interpolation involves evaluating spline basis functions. Grid-based interpolation scales well to large grids because it uses binary search to locate the relevant cell.

Derivatives and Gradients

Grid surfaces support gradient and (for some methods) Hessian evaluation via the Gradient(ReadOnlySpan<Double>, Span<Double>) and Hessian(ReadOnlySpan<Double>, Span2D<Double>) methods.

Derivative behavior depends on the interpolation method:

Nearest-neighbor: Gradient and Hessian are effectively zero (piecewise constant).
Linear: Gradient is constant within each cell, discontinuous at boundaries. Hessian is zero within cells.
Cubic: Gradient and Hessian are continuous and computed analytically from the spline basis functions.
Shape-preserving: Gradient is available; Hessian is not supported for monotone cubic splines.

// Define a 2D grid
double[] xAxis = { 0.0, 1.0, 2.0 };
double[] yAxis = { 0.0, 1.0, 2.0 };

// Array dimension order: values[i, j] = f(xAxis[i], yAxis[j])
double[,] values = new double[xAxis.Length, yAxis.Length];
for (int i = 0; i < xAxis.Length; i++)
{
    for (int j = 0; j < yAxis.Length; j++)
    {
        values[i, j] = i + j;  // Example: f(x,y) = x + y
    }
}

// Create cubic interpolator which supports gradients
var surface = Numerics.NET.Interpolation.GridCubicInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Throw);

// Compute gradient at a point
Vector<double> gradient = Vector.Create(0.0, 0.0);
surface.Gradient(1.0, 1.0, gradient);
Console.WriteLine($"Gradient at (1.0, 1.0): [{gradient[0]:F3}, {gradient[1]:F3}]");

Visual Basic

' Define a 2D grid
Dim xAxis() As Double = { 0.0, 1.0, 2.0 }
Dim yAxis() As Double = { 0.0, 1.0, 2.0 }

' Array dimension order: values(i, j) = f(xAxis(i), yAxis(j))
Dim values(xAxis.Length - 1, yAxis.Length - 1) As Double
For i As Integer = 0 To xAxis.Length - 1
    For j As Integer = 0 To yAxis.Length - 1
        values(i, j) = i + j  ' Example: f(x,y) = x + y
    Next
Next

' Create cubic interpolator which supports gradients
Dim surface = Numerics.NET.Interpolation.GridCubicInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Throw)

' Compute gradient at a point
Dim gradient = Vector.Create(0.0, 0.0)
surface.Gradient(1.0, 1.0, gradient)
Console.WriteLine($"Gradient at (1.0, 1.0): [{gradient(0):F3}, {gradient(1):F3}]")

Visual Basic

No code example is currently available or this language may not be supported.

// Define a 2D grid
let xAxis = [| 0.0; 1.0; 2.0 |]
let yAxis = [| 0.0; 1.0; 2.0 |]

// Array dimension order: values.[i, j] = f(xAxis.[i], yAxis.[j])
let values = Array2D.init xAxis.Length yAxis.Length (fun i j -> 
    float i + float j)  // Example: f(x,y) = x + y

// Create cubic interpolator which supports gradients
let surface = Numerics.NET.Interpolation.GridCubicInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Throw)

// Compute gradient at a point
let gradient = Vector.Create(0.0, 0.0)
surface.Gradient(1.0, 1.0, gradient)
printfn $"Gradient at (1.0, 1.0): [{gradient.[0]:F3}, {gradient.[1]:F3}]"

Boundary Conditions vs. Extrapolation Modes

The grid interpolation system distinguishes between two related but distinct concepts: boundary conditions and extrapolation modes. Understanding this distinction is crucial for correct usage of tensor-product splines.

Boundary conditions (represented by GridBoundaryCondition):

When: Applied during interpolation stencil assembly (construction/evaluation time)
What: Control how to access grid values when the interpolation kernel needs samples at indices outside [0, N-1]. For example, a 4-point cubic kernel near the edge may require i=-1 or i=N.
Where: Only used by tensor-product splines (GridSurface.CreateCubicTensorSpline and GridSurface.CreateTensorSpline). Not applicable to specialized 2D/3D cubic surfaces or linear interpolation.
Options: Clamp (default), Periodic, Reflect, Fill

Extrapolation modes (represented by ExtrapolationMode):

When: Applied during query evaluation
What: Control what happens when a query coordinate falls outside the grid domain [x₀, xₙ₋₁]. Should we throw an exception, clamp to boundaries, wrap periodically, etc.?
Where: Apply to all grid surface types
Options: Throw (default), Clamp, Constant, Periodic, Reflect, Default

Example scenario:

Consider a 1D grid with 10 points and a cubic kernel (width 4). At the left edge (i=0), the cubic stencil needs values at indices {-1, 0, 1, 2}. The boundary condition determines how to map index -1:

Clamp: Use value at index 0
Periodic: Use value at index 9 (wrap around)
Reflect: Use value at index 1 (mirror)
Fill: Use a constant (e.g., 0.0)

Now suppose you query the surface at coordinate x = -5 (outside the grid domain). The extrapolation mode determines the response:

Throw: Throw ArgumentOutOfRangeException
Clamp: Evaluate at x = x₀ (the left boundary)
Periodic: Wrap x into the grid domain and evaluate
Constant: Return NaN (or user-specified fill value)

Key insight: Boundary conditions affect how the interpolant is constructed and evaluated near edges. Extrapolation modes affect whether and how to evaluate queries outside the domain.

Extrapolation Modes in Detail

The ExtrapolationMode parameter controls how the interpolator handles query points outside the grid domain.

Extrapolation modes

Mode	Behavior	Category
Throw	Throws an ArgumentOutOfRangeException if any coordinate is outside the grid bounds. This is the default.	Global control
Clamp	Clamps each coordinate to the nearest boundary value, effectively extending the grid with constant extrapolation. Creates a plateau effect beyond boundaries.	Global control
Constant	Returns a fill value (typically NaN or a user-specified constant) for out-of-bounds queries without evaluating the interpolant. Useful for masking invalid regions.	Global control
Periodic	Wraps coordinates modulo the grid extent, treating the grid as periodic. Ideal for angular data, periodic signals, or toroidal domains.	Functional extension
Reflect	Reflects coordinates at the grid boundaries, creating a symmetric extension. Useful for symmetric functions or mirror boundary conditions.	Functional extension
Default	Uses the surface's default extrapolation mode, derived from the boundary condition. Periodic boundary → Periodic extrapolation, Reflect boundary → Reflect extrapolation, otherwise Throw.	Defers to surface default

Implementation categories:

Global control modes (Throw, Constant, Clamp): Handled by the grid infrastructure before interpolation. The interpolation method never sees out-of-bounds coordinates.
Functional extension modes (Periodic, Reflect): Handled by the interpolation strategy. They transform query coordinates and allow the interpolation method to extend the function mathematically beyond the domain.

The following example demonstrates the different extrapolation modes:

double[] xAxis = { 0.0, 1.0, 2.0 };
double[] yAxis = { 0.0, 1.0 };

// Array dimension order: values[i, j] = f(xAxis[i], yAxis[j])
double[,] values = new double[xAxis.Length, yAxis.Length];
for (int i = 0; i < xAxis.Length; i++)
{
    for (int j = 0; j < yAxis.Length; j++)
    {
        values[i, j] = i + j * 3;
    }
}

// ExtrapolationMode.Throw: Out-of-range queries throw an exception
var throwSurface = Numerics.NET.Interpolation.GridLinearInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Throw);

try
{
    double z = throwSurface.Evaluate(3.0, 0.5); // Outside grid
}
catch (ArgumentOutOfRangeException)
{
    Console.WriteLine("Out-of-range query with Throw mode raised exception");
}

// ExtrapolationMode.Clamp: Clamp to boundary
var clampSurface = Numerics.NET.Interpolation.GridLinearInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Clamp);
double z1 = clampSurface.Evaluate(3.0, 0.5);
Console.WriteLine($"Clamped value at (3.0, 0.5): {z1}");

// ExtrapolationMode.Fill: Return NaN
var fillSurface = Numerics.NET.Interpolation.GridLinearInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Fill);
double z2 = fillSurface.Evaluate(3.0, 0.5);
Console.WriteLine($"Fill value at (3.0, 0.5): {z2}");

Visual Basic

Dim xAxis() As Double = { 0.0, 1.0, 2.0 }
Dim yAxis() As Double = { 0.0, 1.0 }

' Array dimension order: values(i, j) = f(xAxis(i), yAxis(j))
Dim values(xAxis.Length - 1, yAxis.Length - 1) As Double
For i As Integer = 0 To xAxis.Length - 1
    For j As Integer = 0 To yAxis.Length - 1
        values(i, j) = i + j * 3
    Next
Next

' ExtrapolationMode.Throw: Out-of-range queries throw an exception
Dim throwSurface = Numerics.NET.Interpolation.GridLinearInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Throw)

Try
    Dim z = throwSurface.Evaluate(3.0, 0.5) ' Outside grid
Catch ex As ArgumentOutOfRangeException
    Console.WriteLine("Out-of-range query with Throw mode raised exception")
End Try

' ExtrapolationMode.Nearest: Clamp to boundary
Dim nearestSurface = Numerics.NET.Interpolation.GridLinearInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Nearest)
Dim z1 = nearestSurface.Evaluate(3.0, 0.5)
Console.WriteLine($"Clamped value at (3.0, 0.5): {z1}")

' ExtrapolationMode.Fill: Return NaN
Dim fillSurface = Numerics.NET.Interpolation.GridLinearInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Fill)
Dim z2 = fillSurface.Evaluate(3.0, 0.5)
Console.WriteLine($"Fill value at (3.0, 0.5): {z2}")

Visual Basic

No code example is currently available or this language may not be supported.

let xAxis = [| 0.0; 1.0; 2.0 |]
let yAxis = [| 0.0; 1.0 |]

// Array dimension order: values.[i, j] = f(xAxis.[i], yAxis.[j])
let values = Array2D.init xAxis.Length yAxis.Length (fun i j -> 
    float i + float j * 3.0)

// ExtrapolationMode.Throw: Out-of-range queries throw an exception
let throwSurface = Numerics.NET.Interpolation.GridLinearInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Throw)

try
    let z = throwSurface.Evaluate(3.0, 0.5) // Outside grid
    ()
with
| :? ArgumentOutOfRangeException ->
    printfn "Out-of-range query with Throw mode raised exception"

// ExtrapolationMode.Nearest: Clamp to boundary
let nearestSurface = Numerics.NET.Interpolation.GridLinearInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Nearest)
let z1 = nearestSurface.Evaluate(3.0, 0.5)
printfn $"Clamped value at (3.0, 0.5): {z1}"

// ExtrapolationMode.Fill: Return NaN
let fillSurface = Numerics.NET.Interpolation.GridLinearInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Fill)
let z2 = fillSurface.Evaluate(3.0, 0.5)
printfn $"Fill value at (3.0, 0.5): {z2}"

Note

Extrapolation applies only to the coordinate space. The interpolator does not perform polynomial or spline extrapolation beyond the grid boundaries; it modifies the query coordinates and then interpolates.

Examples

Example 1: Basic 2D interpolation (linear and cubic)

Suppose you have a function f(x, y) = x² + y² sampled on a coarse 2D grid, and you want to evaluate it at arbitrary points. Create a grid interpolator and compare linear and cubic results:

// Define a 2D rectilinear grid
double[] xAxis = { 0.0, 1.0, 2.0, 3.0 };
double[] yAxis = { 0.0, 1.0, 2.0 };

// Define values on the grid
// Array dimension order: values[i, j] = f(xAxis[i], yAxis[j])
double[,] values = new double[xAxis.Length, yAxis.Length];
for (int i = 0; i < xAxis.Length; i++)
{
    for (int j = 0; j < yAxis.Length; j++)
    {
        values[i, j] = i + j * 4;  // Example function
    }
}

// Create a linear grid interpolator
var surface = Numerics.NET.Interpolation.GridLinearInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Throw);

// Evaluate at a single point
double z1 = surface.Evaluate(1.5, 1.0);
Console.WriteLine($"Linear interpolation at (1.5, 1.0): {z1}");

// Evaluate at multiple points
double z2 = surface.Evaluate(0.5, 0.5);
double z3 = surface.Evaluate(2.5, 1.5);
Console.WriteLine($"Values: ({z2:F2}, {z3:F2})");

Visual Basic

' Define a 2D rectilinear grid
Dim xAxis() As Double = { 0.0, 1.0, 2.0, 3.0 }
Dim yAxis() As Double = { 0.0, 1.0, 2.0 }

' Define values on the grid
' Array dimension order: values(i, j) = f(xAxis(i), yAxis(j))
Dim values(xAxis.Length - 1, yAxis.Length - 1) As Double
For i As Integer = 0 To xAxis.Length - 1
    For j As Integer = 0 To yAxis.Length - 1
        values(i, j) = i + j * 4  ' Example function
    Next
Next

' Create a linear grid interpolator
Dim surface = Numerics.NET.Interpolation.GridLinearInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Throw)

' Evaluate at a single point
Dim z1 = surface.Evaluate(1.5, 1.0)
Console.WriteLine($"Linear interpolation at (1.5, 1.0): {z1}")

' Evaluate at multiple points
Dim z2 = surface.Evaluate(0.5, 0.5)
Dim z3 = surface.Evaluate(2.5, 1.5)
Console.WriteLine($"Values: ({z2:F2}, {z3:F2})")

Visual Basic

No code example is currently available or this language may not be supported.

// Define a 2D rectilinear grid
let xAxis = [| 0.0; 1.0; 2.0; 3.0 |]
let yAxis = [| 0.0; 1.0; 2.0 |]

// Define values on the grid
// Array dimension order: values.[i, j] = f(xAxis.[i], yAxis.[j])
let values = Array2D.init xAxis.Length yAxis.Length (fun i j -> 
    float i + float j * 4.0)  // Example function

// Create a linear grid interpolator
let surface = Numerics.NET.Interpolation.GridLinearInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Throw)

// Evaluate at a single point
let z1 = surface.Evaluate(1.5, 1.0)
printfn $"Linear interpolation at (1.5, 1.0): {z1}"

// Evaluate at multiple points
let z2 = surface.Evaluate(0.5, 0.5)
let z3 = surface.Evaluate(2.5, 1.5)
printfn $"Values: ({z2:F2}, {z3:F2})"

Linear interpolation provides a good approximation with minimal computation. For smoother interpolation with continuous derivatives, use the specialized 2D cubic surface:

// Create a cubic spline interpolator for smoother results
var cubicSurface = Numerics.NET.Interpolation.GridCubicInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Throw);

// Evaluate at the same points - cubic will be smoother
double zc1 = cubicSurface.Evaluate(1.5, 1.0);
double zc2 = cubicSurface.Evaluate(0.5, 0.5);
Console.WriteLine($"Cubic interpolation at (1.5, 1.0): {zc1}");
Console.WriteLine($"Cubic interpolation at (0.5, 0.5): {zc2}");

Visual Basic

' Create a cubic spline interpolator for smoother results
Dim cubicSurface = Numerics.NET.Interpolation.GridCubicInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Throw)

' Evaluate at the same points - cubic will be smoother
Dim zc1 = cubicSurface.Evaluate(1.5, 1.0)
Dim zc2 = cubicSurface.Evaluate(0.5, 0.5)
Console.WriteLine($"Cubic interpolation at (1.5, 1.0): {zc1}")
Console.WriteLine($"Cubic interpolation at (0.5, 0.5): {zc2}")

Visual Basic

No code example is currently available or this language may not be supported.

// Create a cubic spline interpolator for smoother results
let cubicSurface = Numerics.NET.Interpolation.GridCubicInterpolator(
    xAxis, yAxis, values, ExtrapolationMode.Throw)

// Evaluate at the same points - cubic will be smoother
let zc1 = cubicSurface.Evaluate(1.5, 1.0)
let zc2 = cubicSurface.Evaluate(0.5, 0.5)
printfn $"Cubic interpolation at (1.5, 1.0): {zc1}"
printfn $"Cubic interpolation at (0.5, 0.5): {zc2}"

Example 2: Tensor-product cubic splines with custom boundary conditions

For N-dimensional grids with custom boundary conditions, use CreateCubicTensorSpline or CreateTensorSpline. This example shows periodic interpolation for angular data:

// Create a 2D grid for latitude-longitude data
var latitudes = Vector.LinSpace(-90, 90, 37);   // -90 to +90 degrees
var longitudes = Vector.LinSpace(0, 360, 73);   // 0 to 360 degrees (periodic)

// Sample function values (e.g., temperature field)
var values = new double[latitudes.Count * longitudes.Count];
for (int i = 0; i < values.Length; i++)
    values[i] = Math.Sin(i * 0.1);  // Example data

// Create tensor-product cubic spline with periodic boundary in longitude direction
var surface = Numerics.NET.Interpolation.GridCubicTensorSplineInterpolator(
    new[] { latitudes, longitudes },
    values,
    boundaryCondition: GridBoundaryCondition.Periodic,
    extrapolation: ExtrapolationMode.Default  // Will use Periodic extrapolation
);

// Query at longitude = 370 degrees (wraps to 10 degrees)
double result = surface.Evaluate(new[] { 45.0, 370.0 });
Console.WriteLine($"Temperature at (45°N, 370°E): {result:F2}");

Visual Basic

' Create a 2D grid for latitude-longitude data
Dim latitudes = Vector.LinSpace(-90, 90, 37)   ' -90 to +90 degrees
Dim longitudes = Vector.LinSpace(0, 360, 73)   ' 0 to 360 degrees (periodic)

' Sample function values (e.g., temperature field)
Dim values = New Double(latitudes.Count * longitudes.Count - 1) {}
For i As Integer = 0 To values.Length - 1
    values(i) = Math.Sin(i * 0.1)  ' Example data
Next

' Create tensor-product cubic spline with periodic boundary in longitude direction
Dim surface = Numerics.NET.Interpolation.GridCubicTensorSplineInterpolator(
    New Vector(Of Double)() { latitudes, longitudes },
    values,
    boundaryCondition := GridBoundaryCondition.Periodic,
    extrapolation := ExtrapolationMode.Default  ' Will use Periodic extrapolation
)

' Query at longitude = 370 degrees (wraps to 10 degrees)
Dim result = surface.Evaluate(New Double() { 45.0, 370.0 })
Console.WriteLine($"Temperature at (45°N, 370°E): {result:F2}")

Visual Basic

No code example is currently available or this language may not be supported.

// Create a 2D grid for latitude-longitude data
let latitudes = Vector.LinSpace(-90, 90, 37)   // -90 to +90 degrees
let longitudes = Vector.LinSpace(0, 360, 73)   // 0 to 360 degrees (periodic)

// Sample function values (e.g., temperature field)
let values = Array.init (latitudes.Count * longitudes.Count) (fun i -> Math.Sin(float i * 0.1))

// Create tensor-product cubic spline with periodic boundary in longitude direction
let surface = Numerics.NET.Interpolation.GridCubicTensorSplineInterpolator(
    [| latitudes; longitudes |],
    values,
    boundaryCondition = GridBoundaryCondition.Periodic,
    extrapolation = ExtrapolationMode.Default  // Will use Periodic extrapolation
)

// Query at longitude = 370 degrees (wraps to 10 degrees)
let result = surface.Evaluate([| 45.0; 370.0 |])
printfn $"Temperature at (45°N, 370°E): {result:F2}"

Example 3: Using different spline kernels

The CreateTensorSpline method accepts custom spline kernels. You can choose between linear and cubic kernels, or implement your own:

// Create a 3D grid
var xAxis = Vector.LinSpace(0, 10, 11);
var yAxis = Vector.LinSpace(0, 10, 11);
var zAxis = Vector.LinSpace(0, 10, 11);
var values = new double[11 * 11 * 11];
for (int i = 0; i < values.Length; i++)
    values[i] = i * 0.01;  // Example data

// Option 1: Linear tensor-product interpolation (equivalent to trilinear)
var linearSurface = GridSurface.CreateTensorSpline(
    new[] { xAxis, yAxis, zAxis },
    values,
    kernel: SplineKernels.Linear
);

// Option 2: Cubic Hermite tensor-product interpolation
var cubicSurface = GridSurface.CreateTensorSpline(
    new[] { xAxis, yAxis, zAxis },
    values,
    kernel: SplineKernels.CubicHermite
);

// Both support gradient evaluation
var grad = cubicSurface.Gradient(new[] { 5.0, 5.0, 5.0 });
Console.WriteLine($"Gradient at (5, 5, 5): [{grad[0]:F4}, {grad[1]:F4}, {grad[2]:F4}]");

Visual Basic

' Create a 3D grid
Dim xAxis = Vector.LinSpace(0, 10, 11)
Dim yAxis = Vector.LinSpace(0, 10, 11)
Dim zAxis = Vector.LinSpace(0, 10, 11)
Dim values = New Double(11 * 11 * 11 - 1) {}
For i As Integer = 0 To values.Length - 1
    values(i) = i * 0.01  ' Example data
Next

' Option 1: Linear tensor-product interpolation (equivalent to trilinear)
Dim linearSurface = GridSurface.CreateTensorSpline(
    New Vector(Of Double)() { xAxis, yAxis, zAxis },
    values,
    kernel := SplineKernels.Linear
)

' Option 2: Cubic Hermite tensor-product interpolation
Dim cubicSurface = GridSurface.CreateTensorSpline(
    New Vector(Of Double)() { xAxis, yAxis, zAxis },
    values,
    kernel := SplineKernels.CubicHermite
)

' Both support gradient evaluation
Dim grad = cubicSurface.Gradient(New Double() { 5.0, 5.0, 5.0 })
Console.WriteLine($"Gradient at (5, 5, 5): [{grad(0):F4}, {grad(1):F4}, {grad(2):F4}]")

Visual Basic

No code example is currently available or this language may not be supported.

// Create a 3D grid
let xAxis = Vector.LinSpace(0, 10, 11)
let yAxis = Vector.LinSpace(0, 10, 11)
let zAxis = Vector.LinSpace(0, 10, 11)
let values = Array.init (11 * 11 * 11) (fun i -> float i * 0.01)

// Option 1: Linear tensor-product interpolation (equivalent to trilinear)
let linearSurface = GridSurface.CreateTensorSpline(
    [| xAxis; yAxis; zAxis |],
    values,
    kernel = SplineKernels.Linear
)

// Option 2: Cubic Hermite tensor-product interpolation
let cubicSurface = GridSurface.CreateTensorSpline(
    [| xAxis; yAxis; zAxis |],
    values,
    kernel = SplineKernels.CubicHermite
)

// Both support gradient evaluation
let grad = cubicSurface.Gradient([| 5.0; 5.0; 5.0 |])
printfn $"Gradient at (5, 5, 5): [{grad.[0]:F4}, {grad.[1]:F4}, {grad.[2]:F4}]"

Example 4: Specialized 2D/3D cubic vs. tensor-product cubic

For 2D and 3D grids, you can choose between specialized cubic surfaces (higher smoothness) and tensor-product cubic splines (more flexibility):

var xAxis = Vector.LinSpace(0, 10, 21);
var yAxis = Vector.LinSpace(0, 10, 21);
var values = new double[21, 21];
for (int i = 0; i < 21; i++)
    for (int j = 0; j < 21; j++)
        values[i, j] = Math.Sin(i * 0.3) * Math.Cos(j * 0.3);

// Flatten for GridSurface factory methods
var flatValues = new double[21 * 21];
for (int i = 0; i < 21; i++)
    for (int j = 0; j < 21; j++)
        flatValues[i * 21 + j] = values[i, j];

// Option A: Specialized 2D cubic surface (C² smoothness, natural boundaries)
var specializedCubic = GridSurface.CreateCubic(
    new[] { xAxis, yAxis },
    flatValues,
    extrapolation: ExtrapolationMode.Throw
);

// Option B: Tensor-product cubic (C¹ smoothness, custom boundaries)
var tensorCubic = GridSurface.CreateCubicTensorSpline(
    new[] { xAxis, yAxis },
    flatValues,
    boundaryCondition: GridBoundaryCondition.Reflect,  // Mirror at boundaries
    extrapolation: ExtrapolationMode.Reflect
);

// specializedCubic has smoother derivatives, tensorCubic is more flexible
double z1 = specializedCubic.Evaluate(new[] { 5.0, 5.0 });
double z2 = tensorCubic.Evaluate(new[] { 5.0, 5.0 });
Console.WriteLine($"Specialized: {z1:F4}, Tensor-product: {z2:F4}");

Visual Basic

Dim xAxis = Vector.LinSpace(0, 10, 21)
Dim yAxis = Vector.LinSpace(0, 10, 21)
Dim values = New Double(20, 20) {}
For i As Integer = 0 To 20
    For j As Integer = 0 To 20
        values(i, j) = Math.Sin(i * 0.3) * Math.Cos(j * 0.3)
    Next
Next

' Flatten for GridSurface factory methods
Dim flatValues = New Double(21 * 21 - 1) {}
For i As Integer = 0 To 20
    For j As Integer = 0 To 20
        flatValues(i * 21 + j) = values(i, j)
    Next
Next

' Option A: Specialized 2D cubic surface (C² smoothness, natural boundaries)
Dim specializedCubic = GridSurface.CreateCubic(
    New Vector(Of Double)() { xAxis, yAxis },
    flatValues,
    extrapolation := ExtrapolationMode.Throw
)

' Option B: Tensor-product cubic (C¹ smoothness, custom boundaries)
Dim tensorCubic = GridSurface.CreateCubicTensorSpline(
    New Vector(Of Double)() { xAxis, yAxis },
    flatValues,
    boundaryCondition := GridBoundaryCondition.Reflect,  ' Mirror at boundaries
    extrapolation := ExtrapolationMode.Reflect
)

' specializedCubic has smoother derivatives, tensorCubic is more flexible
Dim z1 = specializedCubic.Evaluate(New Double() { 5.0, 5.0 })
Dim z2 = tensorCubic.Evaluate(New Double() { 5.0, 5.0 })
Console.WriteLine($"Specialized: {z1:F4}, Tensor-product: {z2:F4}")

Visual Basic

No code example is currently available or this language may not be supported.

let xAxis = Vector.LinSpace(0, 10, 21)
let yAxis = Vector.LinSpace(0, 10, 21)
let values = Array2D.init 21 21 (fun i j -> Math.Sin(float i * 0.3) * Math.Cos(float j * 0.3))

// Flatten for GridSurface factory methods
let flatValues = Array.init (21 * 21) (fun idx ->
    let i = idx / 21
    let j = idx % 21
    values.[i, j]
)

// Option A: Specialized 2D cubic surface (C² smoothness, natural boundaries)
let specializedCubic = GridSurface.CreateCubic(
    [| xAxis; yAxis |],
    flatValues,
    extrapolation = ExtrapolationMode.Throw
)

// Option B: Tensor-product cubic (C¹ smoothness, custom boundaries)
let tensorCubic = GridSurface.CreateCubicTensorSpline(
    [| xAxis; yAxis |],
    flatValues,
    boundaryCondition = GridBoundaryCondition.Reflect,  // Mirror at boundaries
    extrapolation = ExtrapolationMode.Reflect
)

// specializedCubic has smoother derivatives, tensorCubic is more flexible
let z1 = specializedCubic.Evaluate([| 5.0; 5.0 |])
let z2 = tensorCubic.Evaluate([| 5.0; 5.0 |])
printfn $"Specialized: {z1:F4}, Tensor-product: {z2:F4}"

Example 5: High-dimensional lookup table

Tensor-product interpolation scales to arbitrary dimensions, making it ideal for multi-parameter lookup tables:

// Create a 3D grid interpolator
Vector<double> xAxis = Vector.Create(0.0, 1.0, 2.0);
Vector<double> yAxis = Vector.Create(0.0, 1.0);
Vector<double> zAxis = Vector.Create(0.0, 1.0, 2.0);

// Axes are ordered as [z, y, x] for row-major storage
Vector<double>[] axes = new Vector<double>[] { zAxis, yAxis, xAxis };

// Flat values array in row-major order (z varies slowest, x varies fastest)
Vector<double> values = Vector.Create(
    // z=0, y=0: x values
    1.0, 2.0, 3.0,
    // z=0, y=1: x values
    4.0, 5.0, 6.0,
    // z=1, y=0: x values
    7.0, 8.0, 9.0,
    // z=1, y=1: x values
    10.0, 11.0, 12.0,
    // z=2, y=0: x values
    13.0, 14.0, 15.0,
    // z=2, y=1: x values
    16.0, 17.0, 18.0
);

// Create the ND linear interpolator
var field = Numerics.NET.Interpolation.GridLinearInterpolator(
    axes, values, ExtrapolationMode.Throw);

// Evaluate at a 3D point
Vector<double> point = Vector.Create(1.0, 0.5, 1.5);
double value = field.Evaluate(point);
Console.WriteLine($"3D interpolation at (1.0, 0.5, 1.5): {value}");

Visual Basic

' Create a 3D grid interpolator
Dim xAxis = Vector.Create(0.0, 1.0, 2.0)
Dim yAxis = Vector.Create(0.0, 1.0)
Dim zAxis = Vector.Create(0.0, 1.0, 2.0)

' Axes are ordered as [z, y, x] for row-major storage
Dim axes = New Vector(Of Double)() { zAxis, yAxis, xAxis }

' Flat values array in row-major order (z varies slowest, x varies fastest)
Dim values = Vector.Create(
    1.0, 2.0, 3.0,
    4.0, 5.0, 6.0,
    7.0, 8.0, 9.0,
    10.0, 11.0, 12.0,
    13.0, 14.0, 15.0,
    16.0, 17.0, 18.0
)

' Create the ND linear interpolator
Dim field = Numerics.NET.Interpolation.GridLinearInterpolator(
    axes, values, ExtrapolationMode.Throw)

' Evaluate at a 3D point
Dim point = Vector.Create(1.0, 0.5, 1.5)
Dim value = field.Evaluate(point)
Console.WriteLine($"3D interpolation at (1.0, 0.5, 1.5): {value}")

Visual Basic

No code example is currently available or this language may not be supported.

// Create a 3D grid interpolator
let xAxis = Vector.Create(0.0, 1.0, 2.0)
let yAxis = Vector.Create(0.0, 1.0)
let zAxis = Vector.Create(0.0, 1.0, 2.0)

// Axes are ordered as [z, y, x] for row-major storage
let axes = [| zAxis; yAxis; xAxis |]

// Flat values array in row-major order (z varies slowest, x varies fastest)
let values = Vector.Create(
    1.0, 2.0, 3.0,
    4.0, 5.0, 6.0,
    7.0, 8.0, 9.0,
    10.0, 11.0, 12.0,
    13.0, 14.0, 15.0,
    16.0, 17.0, 18.0
)

// Create the ND linear interpolator
let field = Numerics.NET.Interpolation.GridLinearInterpolator(
    axes, values, ExtrapolationMode.Throw)

// Evaluate at a 3D point
let point = Vector.Create(1.0, 0.5, 1.5)
let value = field.Evaluate(point)
printfn $"3D interpolation at (1.0, 0.5, 1.5): {value}"

This approach works efficiently for moderate dimensions (typically up to 4-6 dimensions in practice, depending on grid resolution and performance requirements).

Interpolation on rectilinear grids

Rectilinear Grids

Supported Interpolation Methods

Specialized 2D/3D Cubic vs. Tensor-Product Splines

Creating Grid-Based Interpolators

Array Dimension Order

2D and 3D Examples

Evaluating Grid Surfaces

Derivatives and Gradients

Boundary Conditions vs. Extrapolation Modes

Extrapolation Modes in Detail

Examples

See Also

Reference

Other Resources