b样条曲线拟合c++

以下是基于C++的三次B样条曲线拟合代码，包含插值拟合和近似拟合：

#include <iostream>
#include <vector>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

// 三次B样条基函数
double B(int i, int k, double t, const vector<double>& knots) {
    if (k == 0) {
        if (t >= knots[i] && t < knots[i+1]) {
            return 1.0;
        } else {
            return 0.0;
        }
    } else {
        double w1 = 0.0, w2 = 0.0;
        if (knots[i+k] != knots[i]) {
            w1 = (t - knots[i]) / (knots[i+k] - knots[i]) * B(i, k-1, t, knots);
        }
        if (knots[i+k+1] != knots[i+1]) {
            w2 = (knots[i+k+1] - t) / (knots[i+k+1] - knots[i+1]) * B(i+1, k-1, t, knots);
        }
        return w1 + w2;
    }
}

// 三次B样条曲线拟合
void cubicBSplineFit(const vector<double>& x, const vector<double>& y, int n, int k, const vector<double>& knots, vector<double>& c) {
    MatrixXd A(n, n);
    VectorXd b(n);
    for (int i = 0; i < n; i++) {
        b(i) = 0.0;
        for (int j = 0; j < k; j++) {
            if (i+j < n) {
                A(i, i+j) = B(i, k-1, knots[i+j], knots);
            }
        }
    }
    for (int i = 0; i < k; i++) {
        A(n-1, i) = B(n-1, k-1, knots[i], knots);
    }
    b(n-1) = y[n-1];
    c = A.colPivHouseholderQr().solve(b);
}

int main() {
    // 插值拟合
    vector<double> x = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0};
    vector<double> y = {0.0, 0.8, 0.9, 0.1, -0.8, -1.0, -0.6, 0.4, 1.0, 0.6};
    int n = x.size();
    int k = 4;
    vector<double> knots(n+k);
    for (int i = 0; i < n+k; i++) {
        if (i < k) {
            knots[i] = x[0];
        } else if (i >= n) {
            knots[i] = x[n-1];
        } else {
            knots[i] = x[i-k+1];
        }
    }
    vector<double> c;
    cubicBSplineFit(x, y, n, k, knots, c);
    cout << "插值拟合结果：" << endl;
    for (int i = 0; i < n-1; i++) {
        for (double t = knots[i]; t < knots[i+1]; t += 0.1) {
            double y = 0.0;
            for (int j = 0; j < k; j++) {
                y += c[i+j] * B(i, k-1, t, knots);
            }
            cout << t << " " << y << endl;
        }
    }

    // 近似拟合
    vector<double> x2 = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0};
    vector<double> y2 = {0.0, 0.8, 0.9, 0.1, -0.8, -1.0, -0.6, 0.4, 1.0, 0.6};
    int n2 = x2.size();
    int k2 = 4;
    vector<double> knots2(n2+k2);
    for (int i = 0; i < n2+k2; i++) {
        if (i < k2) {
            knots2[i] = x2[0];
        } else if (i >= n2) {
            knots2[i] = x2[n2-1];
        } else {
            knots2[i] = x2[i-k2+1];
        }
    }
    vector<double> c2;
    cubicBSplineFit(x2, y2, n2, k2, knots2, c2);
    cout << "近似拟合结果：" << endl;
    for (int i = 0; i < n2-1; i++) {
        for (double t = knots2[i]; t < knots2[i+1]; t += 0.1) {
            double y = 0.0;
            for (int j = 0; j < k2; j++) {
                y += c2[i+j] * B(i, k2-1, t, knots2);
            }
            cout << t << " " << y << endl;
        }
    }

    return 0;
}