Better way than if else if else... for linear interpolation

NoSenseEtAl picture NoSenseEtAl · Jul 9, 2012 · Viewed 15.5k times · Source

question is easy. Lets say you have function

double interpolate (double x);

and you have a table that has map of known x-> y
for example
5 15
7 18
10 22
note: real tables are bigger ofc, this is just example.

so for 8 you would return 18+((8-7)/(10-7))*(22-18)=19.3333333

One cool way I found is http://www.bnikolic.co.uk/blog/cpp-map-interp.html (long story short it uses std::map, key= x, value = y for x->y data pairs).

If somebody asks what is the if else if else way in title it is basically:

if ((x>=5) && (x<=7))
{
//interpolate
}
else 
     if((x>=7) && x<=10)
     {
      //interpolate
     }

So is there a more clever way to do it or map way is the state of the art? :)

Btw I prefer soutions in C++ but obviously any language solution that has 1:1 mapping to C++ is nice.

Answer

Daniel Fleischman picture Daniel Fleischman · Jul 26, 2012

Well, the easiest way I can think of would be using a binary search to find the point where your point lies. Try to avoid maps if you can, as they are very slow in practice.

This is a simple way:

const double INF = 1.e100;
vector<pair<double, double> > table;
double interpolate(double x) {
    // Assumes that "table" is sorted by .first
    // Check if x is out of bound
    if (x > table.back().first) return INF;
    if (x < table[0].first) return -INF;
    vector<pair<double, double> >::iterator it, it2;
    // INFINITY is defined in math.h in the glibc implementation
    it = lower_bound(table.begin(), table.end(), make_pair(x, -INF));
    // Corner case
    if (it == table.begin()) return it->second;
    it2 = it;
    --it2;
    return it2->second + (it->second - it2->second)*(x - it2->first)/(it->first - it2->first);
}
int main() {
    table.push_back(make_pair(5., 15.));
    table.push_back(make_pair(7., 18.));
    table.push_back(make_pair(10., 22.));
    // If you are not sure if table is sorted:
    sort(table.begin(), table.end());
    printf("%f\n", interpolate(8.));
    printf("%f\n", interpolate(10.));
    printf("%f\n", interpolate(10.1));
}