mlinterp

mlinterp is a fast C++ routine for linear interpolation in arbitrary dimensions (i.e., multilinear interpolation).

mlinterp is written by Parsiad Azimzadeh and released under a permissive MIT License. The latest release can be downloaded here.

Installing

Quick and dirty

mlinterp is a header-only library, meaning that the fastest way to use it is to copy the file mlinterp/mlinterp.hpp into your own project directory and include it using #include <mlinterp.hpp> .

Clean install

If you are on a UNIX-like system with CMake installed, you can download mlinterp and run the following commands (from the mlinterp project directory) to install it to your system:

cmake . sudo make install # Run this to install make # Run this to compile examples

Examples

1d

Let’s interpolate y = sin(x) on the interval [-pi, pi] using 15 evenly-spaced data points.

using namespace mlinterp ; // Boundaries of the interval [-pi, pi] constexpr double b = 3.14159265358979323846 , a = - b ; // Subdivide the interval [-pi, pi] using 15 evenly-spaced points and // evaluate sin(x) at each of those points constexpr int nxd = 15 , nd [] = { nxd }; double xd [ nxd ]; double yd [ nxd ]; for ( int n = 0 ; n < nxd ; ++ n ) { xd [ n ] = a + ( b - a ) / ( nxd - 1 ) * n ; yd [ n ] = sin ( xd [ n ]); } // Subdivide the interval [-pi, pi] using 100 evenly-spaced points // (these are the points at which we interpolate) constexpr int ni = 100 ; double xi [ ni ]; for ( int n = 0 ; n < ni ; ++ n ) { xi [ n ] = a + ( b - a ) / ( ni - 1 ) * n ; } // Perform the interpolation double yi [ ni ]; // Result is stored in this buffer interp ( nd , ni , // Number of points yd , yi , // Output axis (y) xd , xi // Input axis (x) ); // Print the interpolated values cout << scientific << setprecision ( 8 ) << showpos ; for ( int n = 0 ; n < ni ; ++ n ) { cout << xi [ n ] << " \t " << yi [ n ] << endl ; }

2d

Let’s interpolate z = sin(x)cos(y) on the interval [-pi, pi] X [-pi, pi] using 15 evenly-spaced points along the x axis and 15 evenly-spaced points along the y axis.

using namespace mlinterp ; // Boundaries of the interval [-pi, pi] constexpr double b = 3.14159265358979323846 , a = - b ; // Discretize the set [-pi, pi] X [-pi, pi] using 15 evenly-spaced // points along the x axis and 15 evenly-spaced points along the y axis // and evaluate sin(x)cos(y) at each of those points constexpr int nxd = 15 , nyd = 15 , nd [] = { nxd , nyd }; double xd [ nxd ]; for ( int i = 0 ; i < nxd ; ++ i ) { xd [ i ] = a + ( b - a ) / ( nxd - 1 ) * i ; } double yd [ nyd ]; for ( int j = 0 ; j < nyd ; ++ j ) { yd [ j ] = a + ( b - a ) / ( nyd - 1 ) * j ; } double zd [ nxd * nyd ]; for ( int i = 0 ; i < nxd ; ++ i ) { for ( int j = 0 ; j < nyd ; ++ j ) { const int n = j + i * nyd ; zd [ n ] = sin ( xd [ i ]) * cos ( yd [ j ]); } } // Subdivide the set [-pi, pi] X [-pi, pi] using 100 evenly-spaced // points along the x axis and 100 evenly-spaced points along the y axis // (these are the points at which we interpolate) constexpr int m = 100 , ni = m * m ; double xi [ ni ]; double yi [ ni ]; for ( int i = 0 ; i < m ; ++ i ) { for ( int j = 0 ; j < m ; ++ j ) { const int n = j + i * m ; xi [ n ] = a + ( b - a ) / ( m - 1 ) * i ; yi [ n ] = a + ( b - a ) / ( m - 1 ) * j ; } } // Perform the interpolation double zi [ ni ]; // Result is stored in this buffer interp ( nd , ni , // Number of points zd , zi , // Output axis (z) xd , xi , yd , yi // Input axes (x and y) ); // Print the interpolated values cout << scientific << setprecision ( 8 ) << showpos ; for ( int n = 0 ; n < ni ; ++ n ) { cout << xi [ n ] << " \t " << yi [ n ] << " \t " << zi [ n ] << endl ; }

Orders

In the 2d example, the interp routine assumes that the array zd is “ordered” in a certain way.

In particular, if i is the index associated with the x axis and j is the index associated with the y axis, n = j + i * nyd gives the index of the corresponding element in the zd array.

This is referred to as the natural order, and it generalizes to arbitrary dimensions. If instead we wanted n = i + j * nxd , we would use the reverse natural order, which can be invoked by using instead

interp < rnatord > ( nd , ni , zd , zi , xd , xi , yd , yi );

Custom orders

You are free to define your own ordering scheme. To get a sense of how to do this, it’s easiest to study the code for the natural order: