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

Question

I am try to use the Eigen library to create splines. However once I create a spline, I don't know how to get what the value would be at a given point x.

See example below with explanations of my intentions:

#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)
}

Answer 1

I see how that could be confusing. The Eigen Spline fitting module, as you use it, does not model a function R -> R. You could, for example, build a spiral with it. This means that you cannot expect to get an Y value from an X value, and instead you pick out points on the spline by how far along the spline they are (hence the chord lengths).

It is possible to use the module to model a function, albeit not terribly intuitive: consider your Y values points in R¹, and instead of letting Eigen calculate the chord lengths, supply your own set of knot parameters that are spaced like your X values (scaled down to [0,1] so the algorithm can cope). It could be packaged something like this:

#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;
}