Интерполяция собственного сплайна. Как получить значение y сплайна в произвольной точке x?

Я пытаюсь использовать библиотеку 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)
}

c++ spline eigen3
person windenergy    schedule 23.04.2015    source источник

Ответы (2)


arrow_upward
13
arrow_downward

Я вижу, как это может сбивать с толку. Модуль подгонки Eigen Spline, используемый вами, не моделирует функцию R -> R. Вы можете, например, построить с его помощью спираль. Это означает, что вы не можете ожидать получить значение Y из значения X, и вместо этого вы выбираете точки на сплайне по тому, как далеко они расположены вдоль сплайна (отсюда и длина хорды).

Можно использовать модуль для моделирования функции, хотя это и не очень интуитивно понятно: рассмотрите точки значений Y в R1 и вместо того, чтобы позволить 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;
}
person Wintermute    schedule 23.04.2015
comment
@Wintermute Спасибо за этот ответ, который я нашел, когда искал поддержку сплайнов Eigen. Однако я не понимаю, почему сплайн-интерполяция Эйгена не дает тех же результатов, что и простой естественный сплайн-код Vanilla в C. Кроме того, я не понимаю, почему необходимо масштабирование значений: я пробовал с немасштабированными значениями, и интерполяция, кажется, работает , хотя все еще отличается от простых числовых рецептов C-кода - person volatile; 19.08.2015
comment
@volatile В модуле Eigen Splines используются B-сплайны (ссылка на документацию). Требуемое масштабирование значений, вероятно, является артефактом этого; Придется копаться в деталях. В любом случае, без него интерполируются фиктивные значения. Например, вы обнаружите, что если вы замените scaled_value в приведенном выше коде на return x;, s(30) станет -496.714, а не 17. - person Wintermute; 24.08.2015
comment
Как я могу использовать это в трехмерном пространстве? - person Soccertrash; 21.11.2017

arrow_upward
2
arrow_downward

Решение @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;

требует ~ 2 ГБ ОЗУ и требует 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 );

Для этого требуется несколько микросекунд и требуется очень небольшой объем оперативной памяти.

person acac    schedule 23.10.2020