Fetching contributors… Cannot retrieve contributors at this time
159 lines (125 sloc) 4.74 KB

mlinterp

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

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:

```struct natord {
template <typename Index, Index Dimension>
static Index mux(const Index *nd, const Index *indices) {
Index index = 0, product = 1, i = Dimension - 1;
while(true) {
index += indices[i] * product;
if(i == 0) { break; }
product *= nd[i--];
}
return index;
}
};```
You can’t perform that action at this time.