7
\$\begingroup\$

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?

\$\endgroup\$
3
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.