14

问题很简单。假设你有功能

double interpolate (double x);

并且您有一个具有已知 x-> y 映射的表,
例如
5 15
7 18
10 22
注意:真实表更大 ofc,这只是示例。

所以对于 8 你会返回 18+((8-7)/(10-7))*(22-18)=19.3333333

我发现的一种很酷的方法是 http://www.bnikolic.co.uk/blog/cpp-map-interp.html (长话短说它使用 std::map, key= x, value = y for x->y数据对)。

如果有人问标题中的 if else if else 方式是什么,基本上是:

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

那么有没有更聪明的方法来做到这一点,或者地图方式是最先进的?:)

顺便说一句,我更喜欢 C++ 中的解决方案,但显然任何具有 1:1 映射到 C++ 的语言解决方案都很好。

4

6 回答 6

32

好吧,我能想到的最简单的方法是使用二进制搜索来找到你的观点所在的点。如果可以的话,尽量避免使用地图,因为它们在实践中非常慢。

这是一个简单的方法:

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));
}
于 2012-07-26T17:57:20.627 回答
4

您可以使用二叉搜索树来存储插值数据。当您有大量 N 个插值点时,这很有用,因为插值可以在 O(log N) 时间内执行。但是,在您的示例中,情况似乎并非如此,RedX 建议的线性搜索更合适。

#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);
}
于 2012-07-10T10:36:24.977 回答
4

存储您的点排序:

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

然后从最大值循环到最小值,一旦你低于一个数字,你就知道它是你想要的。

你想这么说6

// 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
于 2012-07-10T09:44:27.543 回答
3

是的,我想你应该在这些区间和自然数之间的映射中思考。我的意思是,只需标记间隔并使用开关:

switch(I) {

    case Int1: //whatever
        break;

      ...

    default:

}

我不知道,这是我首先想到的。

如果您的数字在相对较小的间隔内,则EDIT Switch 比 if-else 更有效(在进行映射时要考虑这一点)

于 2012-07-09T14:25:35.500 回答
3

如果您的 x 坐标必须是不规则间隔的,则按排序顺序存储 x 坐标,并使用二进制搜索查找最近的坐标,例如使用 Daniel Fleischman 的答案。

但是,如果您的问题允许,请考虑对规则间隔的数据进行预插值。所以

   5 15
   7 18
   10 22

变成

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

然后在运行时,您可以使用以下内容插入 O(1):

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

使用

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

这不一定适用于所有问题——你最终可能会失去准确性(如果没有包含所有 x 样本的漂亮网格),如果它会使你的表变得更大,它可能会产生糟糕的缓存惩罚。但对于您开始对 x 坐标有一定控制权的情况,这是一个不错的选择。

于 2016-08-10T17:48:49.590 回答
1

您如何获得它是相当可读和可以理解的,并且对于“聪明”的解决方案有很多话要说。但是,您可以取消下限检查和笨拙&&,因为序列是有序的:

if (x < 5)
  return 0;
else if (x <= 7)
  // interpolate
else if (x <= 10)
  // interpolate
...
于 2012-07-10T09:57:00.207 回答