I've just seen that a unsupported module for interpolation exists in eigen.
I can include it with:

Code: Select all

 #include <unsupported/Eigen/Splines>


Now I have a VectorXd Xcoordinates, a VectorXd Ycoordinates, both made up of 100 points and describing a 2D curve and I wanna know the linear interpolation value given a X value.

What do i have to define then?

There is a SplineFitting<SplineType>::Interpolate method. You'd probably have to play around with that. I didn't find an example either.

I've found this:
http://eigen.tuxfamily.org/dox-devel/unsupported/structEigen_1_1SplineFitting.html

SplineFitting< > function;
function.interpolate(const PointArrayType &pts, DenseIndex degree))

pts are the points in which the curve is interpolated, i guess.
DenseIndex degree, is 1 for linear, 2 for square and so?
But i can't understand if before that I need to define a Spline object.

Sorry guys, I have been off the radar for quite some time. I took a little snippet from the unit tests I have written and hope that explains things a bit better.

Code: Select all

  typedef Spline<double,2> Spline2d;
  typedef Spline2d::PointType PointType;
  typedef Spline2d::KnotVectorType KnotVectorType;
  typedef Spline2d::ControlPointVectorType ControlPointVectorType;

  ControlPointVectorType points = ControlPointVectorType::Random(2,100);
  const Spline2d spline = SplineFitting<Spline2d>::Interpolate(points,3);

  KnotVectorType chord_lengths; // knot parameters
  Eigen::ChordLengths(points, chord_lengths);

  for (Eigen::DenseIndex i=0; i<points.cols(); ++i)
  {
    PointType pt = spline( chord_lengths(i) );
    PointType ref = points.col(i);
    VERIFY( (pt - ref).matrix().norm() < 1e-14 );
  }


There is no need to define a spline in advance. The PointArrayType must be an Eigen. IIRC, it can be either an Eigen::Matrix or an Eigen::Array. For a 2D spline it should have 2 rows and as many points as you want to fit.

Unfortunately, there is not yet any support for approximate spline fitting but any patches would be welcome.

In case you run into troubles, or find bugs, just let me know.

Regards,
Hauke

Hauke wrote:Sorry guys, I have been off the radar for quite some time. I took a little snippet from the unit tests I have written and hope that explains things a bit better.

Code: Select all

  typedef Spline<double,2> Spline2d;
  typedef Spline2d::PointType PointType;
  typedef Spline2d::KnotVectorType KnotVectorType;
  typedef Spline2d::ControlPointVectorType ControlPointVectorType;

  ControlPointVectorType points = ControlPointVectorType::Random(2,100);
  const Spline2d spline = SplineFitting<Spline2d>::Interpolate(points,3);

  KnotVectorType chord_lengths; // knot parameters
  Eigen::ChordLengths(points, chord_lengths);

  for (Eigen::DenseIndex i=0; i<points.cols(); ++i)
  {
    PointType pt = spline( chord_lengths(i) );
    PointType ref = points.col(i);
    VERIFY( (pt - ref).matrix().norm() < 1e-14 );
  }


There is no need to define a spline in advance. The PointArrayType must be an Eigen. IIRC, it can be either an Eigen::Matrix or an Eigen::Array. For a 2D spline it should have 2 rows and as many points as you want to fit.

Unfortunately, there is not yet any support for approximate spline fitting but any patches would be welcome.

In case you run into troubles, or find bugs, just let me know.

Regards,
Hauke

Hi Hauke,

I have started to use the unsupported module of Splines for interpolating in my code and it really works nice.
I would like to know if there is any plan for supporting interpolation of surfaces (at least for 3D) or the module
will only support curves.

noether wrote:I have started to use the unsupported module of Splines for interpolating in my code and it really works nice.
I would like to know if there is any plan for supporting interpolation of surfaces (at least for 3D) or the module
will only support curves.

It's great to hear that the module is working nicely.

I have not yet any plans to incorporate tensor product splines for surface interpolation. I would support volunteers who would like to implement them but that is all I can promise at the moment.

Regrds,
Hauke

I need help with the next stuff:

I have the class Generic_Airplane which will be use Spline for 2D interpolation.

Code: Select all

class Generic_Airplane: public Airplane
{
    private:
        Eigen::Spline<double, 2> _spline;


But the spline will be "populated" once I read some data from a file (so I can not construct it as the earlier post), this takes place within the constructor.
If I define _spline there, just I can not compile it because:

Code: Select all

In file included from generic_airplane.cc:1:0:
generic_airplane.hh: In constructor ‘Generic_Airplane::Generic_Airplane(std::string)’:
generic_airplane.hh:16:32: error: no matching function for call to ‘Eigen::Spline<double, 2>::Spline()’
generic_airplane.hh:16:32: note: candidates are:
In file included from /usr/include/eigen3/unsupported/Eigen/Splines:43:0,
                 from generic_airplane.hh:6,
                 from generic_airplane.cc:1:
/usr/include/eigen3/unsupported/Eigen/src/Splines/Spline.h:82:5: note: template<int OtherDegree> Eigen::Spline::Spline(const Eigen::Spline<_Scalar, Dimension, OtherDegree>&)
/usr/include/eigen3/unsupported/Eigen/src/Splines/Spline.h:82:5: note:   template argument deduction/substitution failed:
In file included from generic_airplane.cc:1:0:
generic_airplane.hh:16:32: note:   candidate expects 1 argument, 0 provided
In file included from /usr/include/eigen3/unsupported/Eigen/Splines:43:0,
                 from generic_airplane.hh:6,
                 from generic_airplane.cc:1:
/usr/include/eigen3/unsupported/Eigen/src/Splines/Spline.h:75:5: note: template<class OtherVectorType, class OtherArrayType> Eigen::Spline::Spline(const OtherVectorType&, const OtherArrayType&)
/usr/include/eigen3/unsupported/Eigen/src/Splines/Spline.h:75:5: note:   template argument deduction/substitution failed:
In file included from generic_airplane.cc:1:0:
generic_airplane.hh:16:32: note:   candidate expects 2 arguments, 0 provided
In file included from /usr/include/eigen3/unsupported/Eigen/Splines:43:0,
                 from generic_airplane.hh:6,
                 from generic_airplane.cc:1:
/usr/include/eigen3/unsupported/Eigen/src/Splines/Spline.h:50:9: note: Eigen::Spline<double, 2>::Spline(const Eigen::Spline<double, 2>&)
/usr/include/eigen3/unsupported/Eigen/src/Splines/Spline.h:50:9: note:   candidate expects 1 argument, 0 provided
make: *** [generic_airplane.o] Error 1


How should I define it in the class?

Hi noether,

the reason is obviously that I am not offering a default constructor. For splines of dynamic degree (as in your case), the solution is easy and I could also add the required default constructor. The idea would be to create a constant spline

Code: Select all

template <typename SplineType>
inline SplineType constantSpline()
{
  // create a degree 0 spline (constant)
  const SplineType::KnotVectorType knots = (SplineType::KnotVectorType(1,2) << 0.0, 1.0).finished();
  const SplineType::ControlPointVectorType ctrls = SplineType::ControlPointVectorType::Zero(2,1);
  return Spline<double,2>(knots, ctrls);
}

and then you could initialize your spline in the Generic_Airplane constructor's initializer list like this

Code: Select all

Generic_Airplane : _spline(constantSpline< Eigen::Spline<double,2> >()) {}

I cannot leave the control point and knot vector members of the spline uninitialized since then you have an invalid object which actually does not represent a spline. We probably do not want this.

I will try to come up with some useful default constructors and for the time being the function I wrote above can be used as a quick workaround.

Regards,
Hauke

I've added the missing default constructor. For splines of dynamic degree, it creates a constant, zero degree spline. For splines of fixed degree, it creates a constant spline of the requested degree. Default splines return 0 at all sites u.

That should fix your issue.

Regards,
Hauke

Hauke wrote:Hi noether,

the reason is obviously that I am not offering a default constructor. For splines of dynamic degree (as in your case), the solution is easy and I could also add the required default constructor. The idea would be to create a constant spline

Code: Select all

template <typename SplineType>
inline SplineType constantSpline()
{
  // create a degree 0 spline (constant)
  const SplineType::KnotVectorType knots = (SplineType::KnotVectorType(1,2) << 0.0, 1.0).finished();
  const SplineType::ControlPointVectorType ctrls = SplineType::ControlPointVectorType::Zero(2,1);
  return Spline<double,2>(knots, ctrls);
}

and then you could initialize your spline in the Generic_Airplane constructor's initializer list like this

Code: Select all

Generic_Airplane : _spline(constantSpline< Eigen::Spline<double,2> >()) {}

I cannot leave the control point and knot vector members of the spline uninitialized since then you have an invalid object which actually does not represent a spline. We probably do not want this.

I will try to come up with some useful default constructors and for the time being the function I wrote above can be used as a quick workaround.

Regards,
Hauke

Thank you very much Hauke, although finally I have solved following the next way:

std::shared_ptr<Eigen::Spline<double, 2>>

And when I have the spline built, just
std::make_shared<Eigen::Spline<double, 2>> (spline);

noether wrote:Thank you very much Hauke, although finally I have solved following the next way:

std::shared_ptr<Eigen::Spline<double, 2>>

And when I have the spline built, just
std::make_shared<Eigen::Spline<double, 2>> (spline);

That's fine too, though you should consider the performance implications of std::shared_ptr (ref counting, thread safetiness machinery, see also [1]). I would suggest to either use nesting by value, as you did in your initial post or to use std::unique_ptr.

Regards,
Hauke

[1] http://herbsutter.com/2012/06/05/gotw-1 ... culty-710/

Hi, Hauke and all involved.
Might be long after the last update here. I am playing around with the Spline class and have a simple question - how do I get a value of the 2d spline (say Spline2d) at the given x, which is not one of the provided points. I need a corresponding ChordLengths value, but have no idea how to get it. Thanks!

Hi A32167,

Code: Select all

Spline<float,2> spline;
auto y = spline(x);

The documentation can be found here: http://eigen.tuxfamily.org/dox-devel/unsupported/classEigen_1_1Spline.html#a785602b502f082dc3e715ffcb89295b4

You do no need a corresponding chord length value in order to evaluate the spline, you can use values of 'x' corresponding to chord length values but do not necessarily have to. Keep in mind, that x has to be from the interval [0;1].

Regards,
Hauke

Hi, Hauke. Thanks for answering!
I still don't quite get it. I'll give an example, using the post above.
I want to use the Spline module for interpolating my data. Here I make a spline over y=x^2 defined on an array of 5 values : 1,2,3,4,6.
While the spline works perfectly well, I don't understand how to get the value at x=5. Iterating through different values of chord values I can see, that for x=5 the spline gives me 25 as expected. Still it's probably not optimal to interpolate this way.

Code: Select all

   Eigen::VectorXd xvals(5);
    xvals << 1,2,3,4,6;
    Eigen::VectorXd yvals(5);
    yvals = xvals.array()*xvals.array();
    std::cout << xvals << std::endl << "--------" << std::endl;
    std::cout << yvals << std::endl << "--------" << std::endl;

    typedef Spline<double,2> Spline2d;
    typedef Spline2d::PointType PointType;
    typedef Spline2d::ControlPointVectorType ControlPointVectorType;

    ControlPointVectorType points(2,5);
    points.row(0) = xvals;
    points.row(1) = yvals;
    const Spline2d spline = SplineFitting<Spline2d>::Interpolate(points,3);

   for (RealType u=0; u<=1; u+=0.025) {
        PointType y = spline(u);
        std::cout << "(" << u << ":" << y(0,0) << "," << y(1,0) << ") ";
    }
    std::cout << std::endl;


That shows

Code: Select all

1
2
3
4
6
--------
1
4
9
16
36
--------
(0:1,1) (0.025:1.32155,1.82879) (0.05:1.60908,2.66712) (0.075:1.86598,3.51395) (0.1:2.09565,4.36821) (0.125:2.30148,5.22885) (0.15:2.48689,6.0948) (0.175:2.65526,6.96501) (0.2:2.80999,7.83843) (0.225:2.95449,8.71399) (0.25:3.09215,9.59063) (0.275:3.22592,10.4674) (0.3:3.35655,11.3442) (0.325:3.48409,12.2211) (0.35:3.60862,13.0979) (0.375:3.7302,13.9748) (0.4:3.84889,14.8518) (0.425:3.96476,15.7289) (0.45:4.07787,16.6062) (0.475:4.18829,17.4837) (0.5:4.29608,18.3613) (0.525:4.40131,19.2392) (0.55:4.50404,20.1173) (0.575:4.60434,20.9958) (0.6:4.70227,21.8745) (0.625:4.7979,22.7536) (0.65:4.89129,23.6331) (0.675:4.9825,24.5129) (0.7:5.07161,25.3932) (0.725:5.15867,26.2739) (0.75:5.24375,27.1552) (0.775:5.32692,28.0369) (0.8:5.40824,28.9192) (0.825:5.48778,29.802) (0.85:5.56559,30.6854) (0.875:5.64175,31.5695) (0.9:5.71632,32.4542) (0.925:5.78937,33.3395) (0.95:5.86095,34.2256) (0.975:5.93114,35.1124)


Hi again,

maybe this code example gives you an idea.

Code: Select all

#include <vector>
#include <iostream>

#include <Eigen/Core>
#include <unsupported/Eigen/Splines>

using namespace Eigen;

double uvalue(double x, double low, double high)
{
  return (x - low)/(high-low);
}

VectorXd uvalues(VectorXd xvals)
{
  const double low = xvals.minCoeff();
  const double high = xvals.maxCoeff();
  for (int i=0; i<xvals.size(); ++i)
  {
    xvals(i) = uvalue(xvals(i), low, high);
  }
  return xvals;
}

int main(int argc, char* argv[])
{
  typedef Spline<double,1> Spline2d;

  const VectorXd xvals = (VectorXd(5) << 1,2,3,4,6).finished();
  const VectorXd yvals = xvals.array().square();
  const Spline2d spline = SplineFitting<Spline2d>::Interpolate(yvals.transpose(), 3, uvalues(xvals).transpose());

  const double step = 0.1;
  for (double x = 1; x < 6 + 0.5*step; x += step)
  {
    std::cout << "(" << x << "," << spline(uvalue(x, xvals.minCoeff(), xvals.maxCoeff())).transpose() << ")\n";
  }
  std::cout << std::endl;
}

The best way to achieve what you want is to specify the knots which should be used for the interpolation by hand.

Your initial code showed a spline which mapped from R^1 to R^2 while the function you wanted to fit the spline to maps from R^1 to R^1. Your way of fitting the curve is fine though I am not sure whether this was what you really wanted.

I hope this helps a bit more,
Hauke