Eigen样条插值 - 如何在任意点x处获取样条y值?

4

我正在尝试使用Eigen库创建样条曲线。然而,一旦我创建了一个样条曲线,我不知道如何在给定的点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)
}
2个回答

16

我理解这可能会让人感到困惑。您使用的Eigen样条拟合模块不能建模一个从R -> R的函数,例如您可以用它来建立一个螺旋线。这意味着您不能期望从X值中得到一个Y值,而是通过沿着样条曲线的距离来选择曲线上的点(因此称为弦长度)。

然而,您可以使用该模块建立一个函数,尽管不是非常直观:将Y值视为R1中的点,并提供自己的节点参数集,这些节点参数按照您的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 谢谢你提供的答案,我在谷歌搜索 Eigen 的样条曲线支持时找到了它。然而,我不明白为什么 Eigen 的样条插值结果与纯粹的 C 代码中的 Vanilla 自然样条曲线不一致。另外,我也不明白为什么需要进行数值缩放:我尝试使用未缩放的值进行插值,似乎也可以工作,只是与纯粹的 vanilla numerical receipes C 代码有所不同。 - volatile
@volatile Eigen Splines模块使用B-Splines(链接到文档)。所需的值缩放可能是这个的一个副产品;我需要去挖掘细节。无论如何,如果没有它,将会插值出虚假的值。例如,您会发现,如果在上面的代码中用return x;替换scaled_value,则s(30)变为-496.714而不是17 - Wintermute
我该如何在三维空间中使用它? - Soccertrash

2

@Wintermute的解决方案对于小向量效果很好。但是如果向量大小很大,它会非常缓慢,并消耗大量的内存。

 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;

需要大约2GB的内存,并在我的系统上需要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。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接