CODEWITHSUNDEEP
c++interpolation
NoSenseEtAl
NoSenseEtAl

Reputation: 30028

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

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.

Upvotes: 14

Views: 17758

Answers (6)

Kevin Holt
Kevin Holt

Reputation: 814

If your x-coordinates must be irregularly spaced, then store the x-coordinates in sorted order, and use a binary search to find the nearest coordinate, for example using Daniel Fleischman's answer.

However, if your problem permits it, consider pre-interpolating to regularly spaced data. So

   5 15
   7 18
   10 22

becomes

   5 15
   6 16.5
   7 18
   8 19.3333333
   9 20.6666667
   10 22

Then at run-time you can interpolate with O(1) using something like this:

double interp1( double x0, double dx, double* y, int n, double xi )
{
   double f = ( xi - x0 ) / dx;
   if (f<0) return y[0];
   if (f>=(n-1)) return y[n-1];
   int i = (int) f;
   double w = f-(double)i;
   return dy[i]*(1.0-w) + dy[i+1]*w;
}

using

double y[6] = {15,16.5,18,19.3333333, 20.6666667, 22 }
double yi = interp1( 5.0 , 1.0 , y, 5, xi );

This isn't necessarily suitable for every problem -- you could end up losing accuracy (if there's no nice grid that contains all your x-samples), and it could have a bad cache penalty if it would make your table much much bigger. But it's a good option for cases where you have some control over the x-coordinates to begin with.

Upvotes: 3

Daniel Fleischman
Daniel Fleischman

Reputation: 657

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));
}

Upvotes: 32

RedX
RedX

Reputation: 15175

Store your points sorted:

index X    Y
1     1 -> 3
2     3 -> 7
3     10-> 8

Then loop from max to min and as soon as you get below a number you know it the one you want.

You want let's say 6 so:

// pseudo
for i = 3 to 1
  if x[i] <= 6
    // you found your range!
    // interpolate between x[i] and x[i - 1]
    break; // Do not look any further
  end
end

Upvotes: 4

Ambroz Bizjak
Ambroz Bizjak

Reputation: 8095

You can use a binary search tree to store the interpolation data. This is beneficial when you have a large set of N interpolation points, as interpolation can then be performed in O(log N) time. However, in your example, this does not seem to be the case, and the linear search suggested by RedX is more appropriate.

#include <stdio.h>
#include <assert.h>

#include <map>

static double interpolate (double x, const std::map<double, double> &table)
{
    assert(table.size() > 0);

    std::map<double, double>::const_iterator it = table.lower_bound(x);

    if (it == table.end()) {
        return table.rbegin()->second;
    } else {
        if (it == table.begin()) {
            return it->second;
        } else {
            double x2 = it->first;
            double y2 = it->second;
            --it;
            double x1 = it->first;
            double y1 = it->second;
            double p = (x - x1) / (x2 - x1);
            return (1 - p) * y1 + p * y2;
        }
    }
}

int main ()
{
    std::map<double, double> table;
    table.insert(std::pair<double, double>(5, 6));
    table.insert(std::pair<double, double>(8, 4));
    table.insert(std::pair<double, double>(9, 5));

    double y = interpolate(5.1, table);

    printf("%f\n", y);
}

Upvotes: 5

acraig5075
acraig5075

Reputation: 10756

How you've already got it is fairly readable and understandable, and there's a lot to be said for that over a "clever" solution. You can however do away with the lower bounds check and clumsy && because the sequence is ordered:

if (x < 5)
  return 0;
else if (x <= 7)
  // interpolate
else if (x <= 10)
  // interpolate
...

Upvotes: 1

synack
synack

Reputation: 1749

Yes, I guess that you should think in a map between those intervals and the natural nummbers. I mean, just label the intervals and use a switch:

switch(I) {

    case Int1: //whatever
        break;

      ...

    default:

}

I don't know, it's the first thing that I thought of.

EDIT Switch is more efficient than if-else if your numbers are within a relative small interval (that's something to take into account when doing the mapping)

Upvotes: 3

Related Questions