用c++编写拟合样条曲线？

huangyhg

继续，写段样例代码用openmp计算

huangyhg

以下是使用Eigen和OpenMP进行并行计算的示例代码，用于拟合1万个点的三次B样条曲线：
#include <iostream>
#include <vector>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
#include <omp.h>

using namespace Eigen;

int main()
{
    const int n = 10000; // 数据点数
    const int k = 4;     // B样条的阶数

    // 生成随机数据
    MatrixXd points = MatrixXd::Random(n, 3);

    // 创建B样条矩阵
    SparseMatrix<double> M(n, n);
    for (int i = 0; i < n; i++) {
        for (int j = i - k + 1; j <= i; j++) {
            if (j < 0 || j >= n) continue;
            double t = (i - j) / (double)k;
            double b = 0;
            if (t >= 0 && t < 1) b = (t*t*t - 3*t*t + 3*t - 1) / 6.0;
            else if (t >= 1 && t < 2) b = (4 - 6*t*t*t + 3*t*t + 3*t - 3*t*t*t*t) / 6.0;
            else if (t >= 2 && t <= 3) b = (t*t*t - 3*t*t + 3*t - 1) / 6.0;
            M.coeffRef(i, j) = b;
        }
    }

    // 计算B样条系数
    SparseMatrix<double> A = M.transpose() * M;
    SparseLU<SparseMatrix<double>> solver;
    solver.analyzePattern(A);
    solver.factorize(A);
    VectorXd coeffs = solver.solve(M.transpose() * points);

    // 计算B样条曲线
    const int npoints = 100; // 插值点数
    MatrixXd interp_points(npoints, 3);
    #pragma omp parallel for
    for (int i = 0; i < npoints; i++) {
        double t = i / (double)(npoints - 1);
        MatrixXd B(1, k);
        for (int j = 0; j < k; j++) {
            double tt = t - j;
            if (tt >= 0 && tt < 1) B(0, j) = (tt*tt*tt - 3*tt*tt + 3*tt - 1) / 6.0;
            else if (tt >= 1 && tt < 2) B(0, j) = (4 - 6*tt*tt*tt + 3*tt*tt + 3*tt - 3*tt*tt*tt*tt) / 6.0;
            else if (tt >= 2 && tt <= 3) B(0, j) = (tt*tt*tt - 3*tt*tt + 3*tt - 1) / 6.0;
        }
        interp_points.row(i) = B * coeffs;
    }

    // 打印结果
    std::cout << "B-spline coefficients: " << coeffs.transpose() << std::endl;
    std::cout << "Interpolated points: " << std::endl << interp_points << std::endl;

// 保存结果到文件
    std:fstream outfile("interp_points.txt");
    if (outfile.is_open()) {
        for (int i = 0; i < npoints; i++) {
            outfile << interp_points(i, 0) << " " << interp_points(i, 1) << " " << interp_points(i, 2) << std::endl;
        }
        outfile.close();
    }

    return 0;
}
在上面的代码中，我们将计算得到的B样条曲线写入到名为“interp_points.txt”的文件中，每行包含一个点的x、y、z坐标。这样，您可以使用其他工具（如MATLAB、Python等）加载数据并进行进一步的处理、可视化等操作。