0

我在使用具有等距输入的直接查找表进行 2D 插值时缺乏性能而苦苦挣扎。

目标是在表 z(x,y) 中找到 4 个值 (z00,z01,z11,z10)，对应于每个输入 x 和 y 的两个最接近的值，(x0 < x < x1) 和 ( y0 < y < y1)。

例如，以下查找：

                                              x
是的 1 2 3
4 20 23 29
6 35 37 43
8 47 50 55

由以下数组表示：

const f32 lookup {20,35,47,23,37,50,29,43,55}

此外，与数组的定义一起，结构提供以下数据以实现更有效的查找：

 lowest_x = 1;
 lowest_y = 4;
 step_x = 1;
 step_y = 2;
 length_x = 3;
 length_y = 3;

该算法最耗时的部分是找到对应于输入 x 和 y 之前和之后值的交集的索引。

我目前的方法是按如下方式计算它们：

鉴于 x0 和 y0 的索引倍数为：

index_x0 = u32 ((x-lowest_x)/step_x);

index_y0 = u32 ((y-lowest_y)/step_y);

那么 x0,x1,y0 和 y1 是：

x0 = lowest_x + index_x0 * step_x ;
x1 = x0 + step_x ;
y0 = lowest_y + index_y0 * step_y ;
y1 = y0 + step_y ;

查找 z(x,y) 中的 4 个必要值是：

   z00_index = index_x0*length_y0+index_y0;

    z00= lookup[z00_index]
    z01= lookup[z00_index+1]
    z10= lookup[z00_index+length_y0];
    z11= lookup[z00_index+length_y0+1];

然后将二维插值作为沿 y0 和 y1 的两个 x 插值和一个 y 插值执行：

zxy0 = (x1-x)/(x1-x0)*z00 + (x-x0)/(x1-x0)*z10;

zxy1 = (x1-x)/(x1-x0)*z01 + (x-x0)/(x1-x0)*z11;

z = (y1-y)/(y1-y0)*zxy0 + (y-y0)/(y-y0)*zxy1;

关于如何改进这一点的任何建议？

4

1 回答 1

0

我不是嵌入式的，但有一些微优化的机会，特别是减少乘法和除法的数量，并避免再次重新计算相同的东西。

如果您首先将索引计算为实数：

double ix = (x - x.base) / x.step;

您将获得一个“真实索引”，您可以将其转换为整数：

unsigned ix0 = ix;

现在这些之间的差异为您提供插值索引：

double r1 = ix - ix0;    // == (x - x0) / (x1 - x0)
double r0 = 1.0 - r1;    // == (x1 - x) / (x1 - x0)

您现在拥有所需的一切x：查找索引和插值系数。您不需要的是x间隔开始和结束时的实际值。y对 get执行相同操作iy0，r0并且r1您的插值是：

double z = r0 * (s0 * z00 + s1 * z01) + r1 * (s0 * z10 + s1 * z11);

这是每个坐标一个除法，您可以通过预先计算因子1.0 / x.step和将其转换为乘法1.0 / y.step。

这是您的原始代码的完整示例实现和建议的改进。在 PC 上测试，我得到了 30%-40% 的加速，但没有在嵌入式上测试。小心：代码不做任何边界检查！

#include <stdlib.h>
#include <stdio.h>

typedef struct Lookup Lookup;
typedef struct Range Range;

struct Range {
    double base;
    double step;
    unsigned length;
};

struct Lookup {
    Range x, y;
    double *data;
    double xdenom;
    double ydenom;
};

double lookup1(const Lookup *lo, double x, double y)
{
    unsigned index_x0 = (x - lo->x.base) / lo->x.step;
    unsigned index_y0 = (y - lo->y.base) / lo->y.step;

    double x0 = lo->x.base + index_x0 * lo->x.step;
    double x1 = x0 + lo->x.step;
    double y0 = lo->y.base + index_y0 * lo->y.step;
    double y1 = y0 + lo->y.step;

    double *p = lo->data + index_x0 * lo->y.length + index_y0;
    double z00 = p[0];
    double z01 = p[1];
    double z10 = p[lo->y.length];
    double z11 = p[lo->y.length + 1];

    double zxy0 = (x1 - x) / (x1 - x0)*z00 + (x - x0) / (x1 - x0)*z10;
    double zxy1 = (x1 - x) / (x1 - x0)*z01 + (x - x0) / (x1 - x0)*z11;

    double z = (y1 - y) / (y1 - y0)*zxy0 + (y - y0) / (y1 - y0)*zxy1;

    return z;
}

double lookup2(const Lookup *lo, double x, double y)
{
    double ix = (x - lo->x.base) * lo->xdenom;
    double iy = (y - lo->y.base) * lo->ydenom;
    
    unsigned ix0 = ix;
    unsigned iy0 = iy;
    
    double r1 = ix - ix0;
    double r0 = 1.0 - r1;
    double s1 = iy - iy0;
    double s0 = 1.0 - s1;

    double *p = lo->data + (ix0 * lo->y.length + iy0);

    double z00 = p[0];
    double z01 = p[1];
    double z10 = p[lo->y.length];
    double z11 = p[lo->y.length + 1];

    double z = r0 * (s0 * z00 + s1 * z01)
             + r1 * (s0 * z10 + s1 * z11);

    return z;
}

double urand(void)
{
    return rand() / (1.0 + RAND_MAX);
}

enum {
    L = 1 << 24
};

int main(void)
{
    Lookup lo = {
        {10, 0.1, 80},
        {10, 0.1, 80},
    };
    
    unsigned i, j;
    double *p;
    unsigned l = L;
    double sum = 0;
    
    lo.xdenom = 1.0 / lo.x.step;
    lo.ydenom = 1.0 / lo.y.step;
    p = lo.data = malloc(lo.x.length * lo.y.length * sizeof(*lo.data));
    
    for (i = 0; i < lo.x.length; i++) {
        for (j = 0; j < lo.y.length; j++) {
            *p++ = urand();
        }
    }
    
    while (l--) {
        double x = lo.x.base + (lo.x.length * lo.x.step) * urand();
        double y = lo.y.base + (lo.y.length * lo.y.step) * urand();
        double z = lookup2(&lo, x, y);
        
        sum += z;
    }
    
    printf("%g\n", sum / L);
    
    free(lo.data);    

    return 0;
}
于 2021-03-10T21:46:11.430 回答