特征样条插值 - 如何在任意点 x 处获取样条曲线 y 值

Eigen Spline interpolation - How to get spline y value at arbitray point x?

本文关键字：获取 样条曲线 任意点 插值 样条 特征      更新时间：2023-10-16

我尝试使用特征库来创建样条曲线。但是，一旦我创建了一个样条曲线，我不知道如何获取给定点 x 的值。

请参阅以下示例，并解释我的意图：

#include <Eigen/Core>
#include <unsupported/Eigen/Splines>
int main(int argc, char const* argv[])
{
    // points at (0,0) (15,12) and (30,17)
    Eigen::MatrixXd points(2, 3);
    points << 0, 15, 30,
              0, 12, 17;
    typedef Eigen::Spline<double, 2> spline2d;
    spline2d s = Eigen::SplineFitting<spline2d>::Interpolate(points, 2);
    // I now have a spline called s.
    // I want to do something like:
    double x = 12.34;
    double new_y = s(x)[1];  // However this s() function uses a chord value. What is a chord value?
    // Is there a:
    double new_y2 = s.eval(x)
}

我明白这可能会令人困惑。特征样条拟合模块在使用时，不会对函数 R -> R 进行建模。例如，你可以用它建立一个螺旋。这意味着您不能期望从 X 值中获取 Y 值，而是通过样条上的点沿样条的距离（因此是弦长度）来挑选点。

可以使用该模块对函数进行建模，尽管不是很直观：考虑 R 1 中的 Y 值点，而不是让 Eigen 计算弦长，而是提供您自己的一组结参数，这些参数的间距类似于您的 X 值（缩小到 [0,1]，以便算法可以应对）。它可以像这样打包：

#include <Eigen/Core>
#include <unsupported/Eigen/Splines>
#include <iostream>
class SplineFunction {
public:
  SplineFunction(Eigen::VectorXd const &x_vec,
                 Eigen::VectorXd const &y_vec)
    : x_min(x_vec.minCoeff()),
      x_max(x_vec.maxCoeff()),
      // Spline fitting here. X values are scaled down to [0, 1] for this.
      spline_(Eigen::SplineFitting<Eigen::Spline<double, 1>>::Interpolate(
                y_vec.transpose(),
                 // No more than cubic spline, but accept short vectors.
                std::min<int>(x_vec.rows() - 1, 3),
                scaled_values(x_vec)))
  { }
  double operator()(double x) const {
    // x values need to be scaled down in extraction as well.
    return spline_(scaled_value(x))(0);
  }
private:
  // Helpers to scale X values down to [0, 1]
  double scaled_value(double x) const {
    return (x - x_min) / (x_max - x_min);
  }
  Eigen::RowVectorXd scaled_values(Eigen::VectorXd const &x_vec) const {
    return x_vec.unaryExpr([this](double x) { return scaled_value(x); }).transpose();
  }
  double x_min;
  double x_max;
  // Spline of one-dimensional "points."
  Eigen::Spline<double, 1> spline_;
};
int main(int argc, char const* argv[])
{
  Eigen::VectorXd xvals(3);
  Eigen::VectorXd yvals(xvals.rows());
  xvals << 0, 15, 30;
  yvals << 0, 12, 17;
  SplineFunction s(xvals, yvals);
  std::cout << s(12.34) << std::endl;
}

@Wintermute解决方案适用于小向量。但是，如果矢量大小很大，则速度非常慢并且会消耗大量RAM。

 Eigen::VectorXd xvals = Eigen::VectorXd::LinSpaced(10000,0.,1);
 Eigen::VectorXd yvals = xvals.array().sin(); 
        
 SplineFunction s(xvals, yvals);
 std::cout << s(0.50000001) << std::endl;

需要 ~2 GB RAM，在我的系统上需要 17 秒（使用了 16 个线程）

为了进行比较，我使用了 gsl

gsl_interp_accel* accel_ptr = gsl_interp_accel_alloc();
gsl_spline* spline_ptr;
spline_ptr = gsl_spline_alloc(  gsl_interp_cspline, 10000 );
gsl_spline_init( spline_ptr, &xvals[0], &yvals[0], 10000 );
std::cout << gsl_spline_eval( spline_ptr, 0.50000001, accel_ptr ) << std::endl;
gsl_spline_free( spline_ptr );
gsl_interp_accel_free( accel_ptr );

这需要一些微秒，并且需要非常少量的RAM

相关文章：