使用GSL实现二维数组插值

meteorology/data, programming/cpp
目录

当我们使用格点数据绘制填充图、等值线图时，就需要使用插值。 C++ 科学计算库 GSL 提供双线性插值算法，可以满足简单应用的需要。

简介

GSL（GNU Scientific Library）是由 GNU 出品的 C/C++ 数值计算开源库，由一系列数值计算程序组成，包括随机数生成器、线性代数、统计、插值等方面。

二维数组插值

给定一组递增的坐标 x1,...,xm，和一组递增的坐标 y1,...,yn，并给定每个二维网格点 (xi,yi) 对应的数据 z(i,j)。 二维数值插值就是一个连续的函数 z(x,y)，满足 z(xi,yi) = z(i,j)。如下图所示：

已知 x 坐标 x1,...,xm 和 y 坐标 y1,...,yn，对应的 n * m 个点，通过插值求任意 x,y 坐标的数值。

插值算法

GSL二维数组插值支持两种插值算法：

  • 双线性插值
  • 双三次插值

下面简要介绍这两种插值算法：

双线性插值

由最近的四个点进行两次线性插值得到结果，如下图所示。

具体公式请参阅《Bilinear interpolation》

双三次插值

由最近的十六个点插值得到结果，如下图所示。

具体公式请参阅《Bicubic interpolation》

GSL插值

已知数据

x 坐标列表，x 坐标从小到大排列

std::vector<double> x_coords

y 坐标列表，y 坐标从大到小排列

std::vector<double> y_coords

格点数据，存放顺序为先 x 后 y。

std::vector<std::vector<double>> matrix

创建插值

设置选用的插值算法

const gsl_interp2d_type *interp_type_ = gsl_interp2d_bilinear;

申请加速器

gsl_interp_accel *interp_x_accel_ = gsl_interp_accel_alloc();
gsl_interp_accel *interp_y_accel_ = gsl_interp_accel_alloc();

为插值申请空间

size_t x_count = x_coords.size();
size_t y_count = y_corrds.size();
spline_ = gsl_spline2d_alloc(interp_type_, x_count, y_count);

申请插值需要的坐标空间和值空间

double *spline_x_coords_ = (double*)calloc(x_count, sizeof(double));
for(int i=0;i<x_count;i++)
{
    spline_x_coords_[i] = x_coords[i];
}
double *spline_y_coords_ = (double*)calloc(y_count, sizeof(double));
for(int i = y_count-1; i >= 0; i--)
{
    spline_y_coords_[ y_count-1-i] = y_coords[i];
}
double *spline_values_ = (double*)calloc(x_count * y_count, sizeof(double));

设置矩阵值，注意输入数据中 y 轴是按从大到小顺序排列的，需要颠倒顺序。

for(int x = 0; x < x_count; x++)
{
    for(int y = y_count-1 ; y >= 0; y--)
    {
        gsl_spline2d_set(spline_, spline_values_, x, y_count - 1 - y , matrix[y][x]);
    }
}

最后初始化插值对象

gsl_spline2d_init(spline_, spline_x_coords_, spline_y_coords_, spline_values_, x_count, y_count);

获取结果

使用 gsl_spline2d_eval 计算某个点的值

double value = gsl_spline2d_eval(spline_, x, y, interp_x_accel_, interp_y_accel_);

释放空间

释放 gsl 插值需要的对象

if(interp_x_accel_)
    gsl_interp_accel_free(interp_x_accel_);
if(interp_y_accel_)
    gsl_interp_accel_free(interp_y_accel_);
if(spline_)
    gsl_spline2d_free(spline_);

释放数组

if(spline_x_coords_)
    free(spline_x_coords_);
if(spline_y_coords_)
    free(spline_y_coords_);
if(spline_values_)
    free(spline_values_);

与 Qwt 结合

我在项目中使用 gsl 库是为了实现在 Qwt 中绘制等值线图和填充图。 使用 gsl 自定义 QwtRasterData 的派生类，实现对二维格点数据的插值。 代码请参看 gist