Linear Interpolation C++

Question

I have to write a collection of methods for performing linear, bilinear and trilinear interpolation. I have also to write some tests to show that interpolation is exact for polynomials (which should be the case using these interpolation methods).

I ended up writing the following classes as core for my interpolation

Here one-second compile/run github repo (gtest required).

#ifndef LINEARINTERPOLATE_H
#define LINEARINTERPOLATE_H
#include<assert.h>
namespace FL{
typedef unsigned int uint;
template< int DIM, typename T = long double>
struct point{
 T coords [DIM] ;
 T val;
 inline const T coord(const int c) const{
 assert(c >= 0 && c < DIM);
 return coords[c];
 }
};
template <class T>
class LinearInterpolator{
public:
 /*
 y
 ^ p1
 | /
 | /
 | /
 | p
 | /
 | p0
 |
 |
 o-------------------------> x
 */
 /* P: point the lie between a and b
 * a,b boundary of the 1D cuboid, a<=b
 */
 point<1,T>& Linear ( point<1,T>& p , const point<1,T>&a, const point<1,T>&b, int c=0 ){
 T x_d = (p.coord(c)-a.coord(c)) * (1/(b.coord(c) - a.coord(c)));
 p.val = Linear(a.val,b.val,x_d);
 return p;
 }
 /* P: point the lie inside the cuboid defined by the first two values of v
 * v[0] <= v[1]
 */
 point<1,T>& Linear ( point<1,T>& p , const point<1,T>* v , int c=0 ){
 T x_d = (p.coord(c)-v[0].coord(c)) * (1/(v[1].coord(c) - v[0].coord(c)));
 p.val = Linear(v[0].val,v[1].val,x_d);
 return p;
 }
/*----------------BILINEAR------------------------
 y
 ^
 | p2.......p3
 | . .
 | . .
 | . .
 | . .
 | p0.......p1
 |
 |
 o-------------------------> x
 p point that lie in the cuboid defined by the 4 values array v
-------------------------------------------------*/
 point<2,T>& Bilinear ( point<2,T>& p , const point<2,T>*v ){
 T x_d = (p.coord(0)-v[0].coord(0)) * (1/(v[1].coord(0) - v[0].coord(0)));
 T y_d = (p.coord(1)-v[0].coord(1)) * (1/(v[2].coord(1) - v[0].coord(1)));
 p.val = Bilinear(v[0].val,v[1].val,v[2].val,v[3].val,x_d,y_d);
 return p;
 }
/*----------------TRILINEAR------------------------
 *
 y p6____________p7
 ^ / | /|
 | / | / |
 | / | / |
 | p2___|_________p3 |
 | | | | |
 | | p4_________|__p5
 | | / | /
 | | / | /
 | | / |/
 | p0_____________p1
 |
 |-------------------------> x
 /
 /
 /
 /
 /
 z/
 p point that lie in the cuboid defined by the 8 values array v
-------------------------------------------------*/
 point<3,T>& Trilinear ( point<3,T>& p , const point<3,T>*v ){
 T x_d = (p.coord(0)-v[0].coord(0)) * (1/(v[1].coord(0) - v[0].coord(0)));
 T y_d = (p.coord(1)-v[0].coord(1)) * (1/(v[2].coord(1) - v[0].coord(1)));
 T z_d = (p.coord(2)-v[0].coord(2)) * (1/(v[4].coord(2) - v[0].coord(2)));
 p.val =Trilinear(v[0].val,v[1].val,v[2].val,v[3].val,v[4].val,v[5].val,v[6].val,v[7].val,
 x_d,y_d,z_d);
 return p;
 }
private:
 inline T Linear(const T f0, const T f1, const T xd){
 //return std::fma(f0, (1-xd) , f1*xd);
 return f0*(1.0-xd) + f1*xd;
 }
 inline T Bilinear(const T f00,const T f10,const T f01, const T f11, const T xd, const T yd){
 const T c0 = f00*(static_cast<T>(1.0)-xd) + f10*xd;
 const T c1 = f01*(static_cast<T>(1.0)-xd) + f11*xd;
 return Linear(c0,c1,yd);
 }
 inline T Trilinear(const T f000, const T f100,const T f010,const T f110,const T f001, const T f101,const T f011,const T f111,const T xd,const T yd ,const T zd){
 const T c0 = Bilinear(f000,f100,f010,f110,xd,yd);
 const T c1 = Bilinear(f001,f101,f011,f111,xd,yd);
 return Linear(c0,c1,zd);
 }
}; //class interpolator
}
#endif /* INTERPOLATE */

Everything is working fine and tests shows that the interpolation is correct for polynomials.

1. Could I improve the design of the class someway?

2. Is there any better way to do the calculation in order to minimize the rounding errors?

Answer 1

I'd prefer to include the C++ header <cassert>, to be consistent. Use the C compatibility headers only in code that must compile with a C compiler.

I don't see the need for the final argument c to Linear(), given that it's passed to point<1,T>::coords(), for which the only valid value is 0. Also there's an inaccurate comment there: a,b boundary of the 1D cuboid, a<=b - but we need a<b to avoid division by zero.

I don't see a good reason for creating a reciprocal and multiplying by that, instead of simply dividing. It appears to me to be unnecessary duplication:

 T x_d = (p.coord(c)-a.coord(c)) * (1/(b.coord(c) - a.coord(c)));

becomes:

 T x_d = (p.coord(c)-a.coord(c)) / (b.coord(c) - a.coord(c));

We might want to make a helper function for these very similar "portion" computations.

We can make point::coord() constexpr, and all the members of LinearInterpolator can be made both constexpr and static.

Instead of the static_cast of 1.0 to T, we might make a simple constant:

private:
 static constexpr T unity = 1.0;

BTW, a "1D cuboid" is normally called a line, and a "2D cuboid" a rectangle.