Using GSL to create a matrix data structure using in Qwt

matrix_data.cpp

#include "matrix_data.h"

MatrixData::MatrixData(QVector<double> x_coords, QVector<double> y_coords, QVector<QVector<double>> matrix):
    QwtRasterData{},
    interp_type_{gsl_interp2d_bilinear},
    interp_x_accel_{nullptr},
    interp_y_accel_{nullptr},
    spline_x_coords_{nullptr},
    spline_y_coords_{nullptr},
    spline_values_{nullptr},
    spline_{nullptr}
{
    interp_x_accel_ = gsl_interp_accel_alloc();
    interp_y_accel_ = gsl_interp_accel_alloc();

    size_t x_count = x_coords.length();
    size_t y_count = y_coords.length();

    spline_ = gsl_spline2d_alloc(interp_type_, x_count, y_count);

    spline_x_coords_ = (double*)calloc(x_count, sizeof(double));
    for(int i=0;i<x_count;i++)
    {
        spline_x_coords_[i] = x_coords[i];
    }

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

    spline_values_ = (double*)calloc(x_count * y_count, sizeof(double));

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

    double max_value = *std::max_element(spline_values_, spline_values_ + x_count*y_count);
    double min_value = *std::min_element(spline_values_, spline_values_ + x_count*y_count);

    setInterval(Qt::XAxis, QwtInterval(x_coords.first(), x_coords.last()));
    setInterval(Qt::YAxis, QwtInterval(y_coords.last(), y_coords.first()));
    setInterval(Qt::ZAxis, QwtInterval(min_value, max_value));
}

MatrixData::~MatrixData()
{
    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_);
}

double MatrixData::value(double x, double y) const
{
    if(spline_ == nullptr)
        return 0;
    double value = gsl_spline2d_eval(spline_, x, y, interp_x_accel_, interp_y_accel_);
    return value;
}

matrix_data.h

#pragma once

#include <qwt/qwt_raster_data.h>
#include <QVector>

#include <gsl/gsl_math.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>

class MatrixData : public QwtRasterData
{
public:
    MatrixData(
            QVector<double> x_coords,
            QVector<double> y_coords,
            QVector<QVector<double>> matrix
    );
    ~MatrixData();

    double value(double x, double y) const override;

private:
    const gsl_interp2d_type *interp_type_;
    gsl_interp_accel *interp_x_accel_;
    gsl_interp_accel *interp_y_accel_;
    gsl_spline2d *spline_;

    double *spline_x_coords_;
    double *spline_y_coords_;
    double *spline_values_;
};