2022.08.08 Monday @BJ

Project 中要对二维数据做个插值，用 C 写了个双线性插值的代码，备忘一下。

// bdr 是旧网格的边界, bdrq 是新网格的边界
// NR 是旧网格的节点数，NRq 是新网格的节点数
// v 是旧数据，vq 是再新网格上的插值结果
void Interp2(double bdr[][2], double *v, int *NR, double bdrq[][2], double *vq, int *NRq)
{
    int DIM = 2;
    int NRD = prod(NR, DIM), NRDq = prod(NRq, DIM);
    double dx = (bdr[0][1] - bdr[0][0]) / NR[0], dy = (bdr[1][1] - bdr[1][0]) / NR[1];
    double dxq = (bdrq[0][1] - bdrq[0][0]) / NRq[0], dyq = (bdrq[1][1] - bdrq[1][0]) / NRq[1];

    double xq, yq;
    int num, index[2], indexq[2];
    double index_x, index_y;
    bool out_of_domain;
    double a, b;
    for (int i = 0; i < NRDq; i++)
    {
        // 计算第 i 个点的坐标
        number2index_C(indexq, i, NRq, DIM);
        xq = bdrq[0][0] + indexq[0] * dxq;
        yq = bdrq[1][0] + indexq[1] * dyq;
        // 计算该点在旧网格中左下角的坐标点
        index_x = (xq - bdr[0][0]) / dx;
        index_y = (yq - bdr[1][0]) / dy;
        index[0] = (int)index_x;
        index[1] = (int)index_y;
        a = index_x - index[0];
        b = index_y - index[1];

        out_of_domain = false;
        if (index_x < 0)
        {
            index[0] = 0;
            out_of_domain = true;
        }
        if (index_y < 0)
        {
            index[1] = 0;
            out_of_domain = true;
        }
        if (index_x >= NR[0] - 1)
        {
            index[0] = NR[0] - 1;
            out_of_domain = true;
        }
        if (index_y >= NR[1] - 1)
        {
            index[1] = NR[1] - 1;
            out_of_domain = true;
        }

        // 超出边界的部分, 用最简单的最近邻插值
        if (out_of_domain)
        {
            num = index2number_C(index, NR, DIM);
            vq[i] = v[num];
            continue;
        }

        // 其余用双线性插值
        vq[i] = (1 - a) * (1 - b) * v[index2number_C(index, NR, DIM)]; // index 为左下角点
        index[0] += 1;
        vq[i] += a * (1 - b) * v[index2number_C(index, NR, DIM)]; // index 为右下角点
        index[1] += 1;
        vq[i] += a * b * v[index2number_C(index, NR, DIM)]; // index 为右上角点
        index[0] -= 1;
        vq[i] += (1 - a) * b * v[index2number_C(index, NR, DIM)]; // index 为左上角点
    }
}

其中用到的指标转化函数为

int index2number_C(int *index, int *N, int D)
{
    int num = index[0];
    for (int i = 1; i<D; i++)
        num = num*N[i] + index[i];
    return num;
}
void number2index_C(int *index, int num, int *N, int D)
{
    for (int i = D - 1; i >= 0; i--)
    {
        index[i] = num % N[i];
        num = (num - index[i]) / N[i];
    }
}