使用GSL实现二维数组插值

目录

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

简介

GSLGNU 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