Uniform Grid Interpolation in 1D, 2D and 3D

Recently I had need for code to interpolate data stored on regular grids in 1D, 2D and 3D.  This code is straightforward but time-consuming to write and I often find myself re-implementing it for each project that I work on.

The linked header file includes templated interpolation routines for 1D, 2D and 3D interpolation on uniformly gridded data with linear, cosine and cubic (Catmull-Rom) interpolants.  It is fully documented with DOXYGEN and provides a simple fairly consistent API.  The code works by interpolating values that are passed in either as individual values or as arrays allowing it to work with arbitrarily defined grids and also optionally returns the weights associated with each value.  I have used the code to implement

The following code snippet taken from one of the 2D test functions shows the use of the API. In all API functions, points are ordered lexographically by axis, i.e. p[0] = [x, y], p[1] = [x+1, y], p[2] = [x, y+1], p[3] = [x+1, y+1].

/**
 @brief test function for the 2D linear interpolation routines, currently just
 checks that the methods are consistent, i.e. that they all produce the same
 output for the same input
 @return true if test passes false otherwise
*/
static bool interpolation_test_lerp_2D(){
 const double t0=0.325f, t1=0.675, v[] = { 1.0, 2.0, 2.0, 5.0 };
 double o0, o1, o2, o3, w[4];
 bool result = true;
 
 // check that API functions are consistent
 o0 = interpolation_lerp_2D( v[0], v[1], v[2], v[3], t0, t1 );
 o1 = interpolation_lerp_2D( v[0], v[1], v[2], v[3], t0, t1, w[0], w[1], w[2], w[3] );
 o2 = interpolation_lerp_2D( v, t0, t1 );
 o3 = interpolation_lerp_2D( v, t0, t1, w );
 result &= fabs(o0-o1) < 1e-9 && fabs(o0-o2) < 1e-9 && fabs(o0-o3) < 1e-9;
 
 // check that input values are reproduced
 result &= fabs( interpolation_lerp_2D( v[0], v[1], v[2], v[3], 0.0, 0.0 )-v[0] ) < 1e-9; 
 result &= fabs( interpolation_lerp_2D( v[0], v[1], v[2], v[3], 1.0, 0.0 )-v[1] ) < 1e-9;
 result &= fabs( interpolation_lerp_2D( v[0], v[1], v[2], v[3], 0.0, 1.0 )-v[2] ) < 1e-9;
 result &= fabs( interpolation_lerp_2D( v[0], v[1], v[2], v[3], 1.0, 1.0 )-v[3] ) < 1e-9; 

 return result;
}

The code is freely available for all purposes and can be downloaded below:

https://sites.google.com/site/jamesgregson/tmp/interpolation.h

Please contact me if you find bugs, questions or find the code useful.