求与Scipy SmoothBivariateSpline等效的C++函数及代码转换方案

原Python代码实现了对非结构化一维numpy数组x、y、z的双变量平滑样条拟合（kx=1, ky=1即双线性），并支持在任意点(mc, yc)求值（含双线性外推）：

dispX = SmoothBivariateSpline(x, y, z, bbox=[0, 2859, 0, 2399],kx=1,ky=1)
result = dispX.ev(mc,yc)    

已尝试Boost（无外推能力）、GSL（仅支持结构化数据）、OpenCV/IPP（无对应样条函数），以下是可行的C++实现方案：

方案一：手动实现Delaunay三角剖分+线性插值/外推 #

kx=1、ky=1的SmoothBivariateSpline本质是基于线性模型的拟合，针对非结构化数据，Delaunay三角剖分是最直接的处理方式，同时可以轻松实现外推：

步骤说明 #

1. 对(x,y)点集做Delaunay三角剖分，将平面划分为多个三角形单元
2. 为每个三角形预计算线性拟合系数（z = a*x + b*y + c）
3. 查询点时：
   • 若点在凸包内，找到所在三角形并代入系数计算
   • 若点在凸包外（需外推），取最近顶点所在的三角形系数进行线性延伸计算

C++代码示例（基于CGAL库） #

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <vector>
#include <iterator>

// 定义CGAL基础类型
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Delaunay_triangulation_2<K> Delaunay;
typedef K::Point_2 Point;

// 存储每个三角形的线性拟合参数
struct TriangleModel {
    double a, b, c; // z = a*x + b*y + c
};

// 预处理：生成三角剖分并计算每个三角形的拟合参数
std::pair<Delaunay, std::vector<TriangleModel>> preprocess(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& z) {
    std::vector<Point> points;
    points.reserve(x.size());
    for (size_t i = 0; i < x.size(); ++i) {
        points.emplace_back(x[i], y[i]);
    }

    Delaunay dt;
    dt.insert(points.begin(), points.end());

    std::vector<TriangleModel> models;
    models.reserve(dt.number_of_finite_faces());

    for (auto face_it = dt.finite_faces_begin(); face_it != dt.finite_faces_end(); ++face_it) {
        // 获取三个顶点的索引
        auto v0 = face_it->vertex(0);
        auto v1 = face_it->vertex(1);
        auto v2 = face_it->vertex(2);
        const size_t idx0 = std::distance(dt.finite_vertices_begin(), v0);
        const size_t idx1 = std::distance(dt.finite_vertices_begin(), v1);
        const size_t idx2 = std::distance(dt.finite_vertices_begin(), v2);

        // 解线性方程组得到拟合参数
        const double x0 = x[idx0], y0 = y[idx0], z0 = z[idx0];
        const double x1 = x[idx1], y1 = y[idx1], z1 = z[idx1];
        const double x2 = x[idx2], y2 = y[idx2], z2 = z[idx2];

        const double det = (x1 - x0) * (y2 - y0) - (y1 - y0) * (x2 - x0);
        const double a = ((z1 - z0) * (y2 - y0) - (z2 - z0) * (y1 - y0)) / det;
        const double b = ((x2 - x0) * (z1 - z0) - (x1 - x0) * (z2 - z0)) / det;
        const double c = z0 - a * x0 - b * y0;

        models.push_back({a, b, c});
    }

    return {std::move(dt), std::move(models)};
}

// 求值函数：支持内插和外推
double evaluate(double mc, double yc, const Delaunay& dt, const std::vector<TriangleModel>& models) {
    const Point query(mc, yc);
    auto face = dt.locate(query);

    if (dt.is_infinite(face)) {
        // 外推：取最近顶点所在的有限三角形进行计算
        auto nearest_vertex = dt.nearest_vertex(query);
        auto face_it = dt.incident_faces(nearest_vertex);
        while (dt.is_infinite(*face_it)) {
            ++face_it;
        }
        const size_t model_idx = std::distance(dt.finite_faces_begin(), face_it);
        const auto& model = models[model_idx];
        return model.a * mc + model.b * yc + model.c;
    } else {
        // 内插：使用所在三角形的拟合模型
        const size_t model_idx = std::distance(dt.finite_faces_begin(), face);
        const auto& model = models[model_idx];
        return model.a * mc + model.b * yc + model.c;
    }
}

方案二：嵌入Python解释器直接调用Scipy #

如果不想手动实现算法，可以在C++中嵌入Python解释器，直接调用原Scipy逻辑，快速完成适配：

C++代码示例 #

#include <Python.h>
#include <numpy/arrayobject.h>
#include <vector>

double scipy_spline_eval(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& z, double mc, double yc) {
    Py_Initialize();
    import_array(); // 初始化numpy模块

    // 将C++数组转换为numpy数组
    const npy_intp dims = static_cast<npy_intp>(x.size());
    PyObject* x_np = PyArray_SimpleNewFromData(1, &dims, NPY_DOUBLE, const_cast<double*>(x.data()));
    PyObject* y_np = PyArray_SimpleNewFromData(1, &dims, NPY_DOUBLE, const_cast<double*>(y.data()));
    PyObject* z_np = PyArray_SimpleNewFromData(1, &dims, NPY_DOUBLE, const_cast<double*>(z.data()));

    // 导入scipy.interpolate模块
    PyObject* interpolate_mod = PyImport_ImportModule("scipy.interpolate");
    PyObject* spline_func = PyObject_GetAttrString(interpolate_mod, "SmoothBivariateSpline");

    // 构造函数参数
    PyObject* args = PyTuple_New(3);
    PyTuple_SetItem(args, 0, x_np);
    PyTuple_SetItem(args, 1, y_np);
    PyTuple_SetItem(args, 2, z_np);

    // 构造关键字参数
    PyObject* kwargs = PyDict_New();
    PyObject* bbox = PyList_New(4);
    PyList_SetItem(bbox, 0, PyLong_FromLong(0));
    PyList_SetItem(bbox, 1, PyLong_FromLong(2859));
    PyList_SetItem(bbox, 2, PyLong_FromLong(0));
    PyList_SetItem(bbox, 3, PyLong_FromLong(2399));
    PyDict_SetItemString(kwargs, "bbox", bbox);
    PyDict_SetItemString(kwargs, "kx", PyLong_FromLong(1));
    PyDict_SetItemString(kwargs, "ky", PyLong_FromLong(1));

    // 创建样条对象并调用ev方法
    PyObject* spline_obj = PyObject_Call(spline_func, args, kwargs);
    PyObject* ev_method = PyObject_GetAttrString(spline_obj, "ev");
    PyObject* ev_args = PyTuple_New(2);
    PyTuple_SetItem(ev_args, 0, PyFloat_FromDouble(mc));
    PyTuple_SetItem(ev_args, 1, PyFloat_FromDouble(yc));
    PyObject* result_obj = PyObject_CallObject(ev_method, ev_args);

    // 转换结果为C++类型
    const double result = PyFloat_AsDouble(result_obj);

    // 释放Python对象资源
    Py_DECREF(result_obj);
    Py_DECREF(ev_args);
    Py_DECREF(ev_method);
    Py_DECREF(spline_obj);
    Py_DECREF(kwargs);
    Py_DECREF(bbox);
    Py_DECREF(args);
    Py_DECREF(spline_func);
    Py_DECREF(interpolate_mod);
    Py_DECREF(x_np);
    Py_DECREF(y_np);
    Py_DECREF(z_np);

    Py_Finalize();
    return result;
}

注意事项 #

• 需要系统安装Python和Scipy库
• 编译时需要链接Python和numpy的动态库（例如：-lpython3.x -lnumpy）

内容的提问来源于stack exchange，提问作者Celeste