ttk592 / spline

c++ cubic spline library

C++ cubic spline interpolation

This is a lightweight implementation of cubic splines to interpolate points f(x_i) = y_i with the following features.

• available spline types:
  • cubic C² splines: global, twice continuously differentiable
  • cubic Hermite splines: local, continuously differentiable (C¹)
• boundary conditions: first and second order derivatives can be specified, not-a-knot condition, periodic condition is not implemented
• extrapolation
  • linear: if first order derivatives are specified or 2nd order = 0
  • quadratic: if 2nd order derivatives not equal to zero specified
• monotonicity can be enforced (when input is monotonic as well)

Usage

The library is a header-only file with no external dependencies and can be used like this:

#include <vector>
#include "spline.h"
...
    std::vector<double> X, Y;
    ...
    // default cubic spline (C^2) with natural boundary conditions (f''=0)
    tk::spline s(X,Y);			// X needs to be strictly increasing
    double value=s(1.3);		// interpolated value at 1.3
    double deriv=s.deriv(1,1.3);	// 1st order derivative at 1.3
    std::vector<double> solutions = s.solve(0.0);	// solves s(x)=0.0
    ...

The constructor can take more arguments to define the spline, e.g.:

    // cubic Hermite splines (C^1) with enforced monotonicity and
    // left curvature equal to 0.0 and right slope equal 1.0
    tk::spline s(X,Y,tk::spline::cspline_hermite, true,
                 tk::spline::second_deriv, 0.0,
                 tk::spline::first_deriv, 1.0);

This is identical to (must be called in that order):

    tk::spline s;
    s.set_boundary(tk::spline::second_deriv, 0.0,
                   tk::spline::first_deriv, 1.0);
    s.set_points(X,Y);
    s.make_monotonic();

Spline types

Splines are piecewise polynomial functions to interpolate points (x_i, y_i). In particular, cubic splines can be represented as

• f(x) = a_i + b_i (x-x_i) + c_i (x-x_i)² + d_i (x-x_i)³, for all x in [x_i, x_i+1)
• f(x_i)=y_i

The following splines are available.

• tk::spline::cspline: cubic C² spline
  • twice continuously differentiable, e.g. f'(x_i) and f''(x_i) exist
  • this, together with boundary conditions uniquely determines the spline
  • requires solving a sparse equation system
  • is a global spline in the sense that changing an input point will impact the spline everywhere
  • setting first order derivatives at the boundary will break C² at the boundary
• tk::spline::cspline_hermite: cubic Hermite spline
  • once continuously differentiable (C¹)
  • first order derivatives are specified by finite differences, e.g. on a uniform x-grid:
    • f'(x_i) = (y_i+1-y_i-1)/(x_i+1-x_i-1)
  • is a local spline in the sense that changing grid points will only affect the spline around that grid point in a few adjacent segments

A function to enforce monotonicity is available as well:

• tk::spline::make_monotonic(): will make the spline monotonic if input grid points are monotonic
  • this function can only be called after set_points(...) has been called
  • it will break C² if the original spline was C² and not already monotonic
  • it will break boundary conditions if it was not monotonic in the first or last segment

References

https://kluge.in-chemnitz.de/opensource/spline/

spline's Issues

How to obtain the curve length of a spline

I need to find the curve length , could you please give me some suggestoins.

x need to increase monotonically

I use this to smooth the path, but the x of the path point is not necessarily monotonically increasing. How should I deal with it?

[question] Interpolating closed spline ?

Hello,

Is this library able to interpolate a closed spline curve ?

Thanks

How to do a spline and divide the spline into points of equal distance?

So task like the following:
4 points input
100 points output in between the first and the last one

first derivative

Is anyone adding functionality to calculate the first derivative at the interior points?

Computation time

Hello. Thank you for this library. It is quite easy to use.
Your benchmarks shows that it takes 0.2sec to access spline 10^7 times.
But it takes 100 seconds in my case. I build a spline and include it as a member to my class.
Then in another method I access this spline about 10^7 times.
Please suggest an algorithm optimization.

More liberal license?

May I suggest a change to a more liberal license (MIT, LGPL, Apache2...)? I see it more appropriate for the scope and the amount of code. Of course it's your choice.

Is this should be h^2 for first_deriv?

spline/src/spline.h

Line 539 in 5894bea

interpol=(3.0*m_d[idx]*h + 2.0*m_c[idx])*h + m_b[idx];

Is this should be h^2 for first_deriv?
interpol = (3.0 * m_d[idx] * h * h + 2.0 * m_c[idx]) * h + m_b[idx];

warning: ‘CLASS’ has a field ‘CLASS::spline_’ whose type uses the anonymous namespace

Usage of spline as a member of class creates warning.

#include "spline/spline.h"
class CLASS {
tk::spline spline_;
}

Can implement a 'not-a-knot' boundary condition like Python?

I found that the default result is the same as Python(scipy.interpolate.CubicSpline)'s natural boundary condition.
But I don't know how to match the default result of Python(with default boundary condition not-a-knot)

Compile error

When running make in the examples folder I get:

g++ -I../../alglib/src -I../src -O2 -DNDEBUG -Wall -c interpol_alglib.cpp -o interpol_alglib.o
interpol_alglib.cpp:4:37: fatal error: interpolation.h: No such file or directory
 #include "interpolation.h" // alglib

Looks like there is a missing file from alglib - perhaps document this better?

Is it a Cubic B-Spline Function?

I ran the code twice with 5 control points, in the second run I changed the fourth control point by a certain magnitude and observed that the global shape of the spline changed completely. I was hoping for more local control.

Not-a-knot index error

Hi @ttk592,

I think this line is wrong: https://github.com/ttk592/spline/blob/master/src/spline.h#L332

Shouldn't it be: rhs[n-1] = 0.0;?

.Solve(-) not working, gives out empty vector

Hi, great work on the splines code. I just have a small issue i stumbled upon using it in my project. I'm trying to recreate this splines into my code. Thanks in advance! Here is the issue, the splines and my code:

 int WorldGen(float Cont, float Eros, float PV ) {
			int worldHeight = 0;
			//Peakes & Valleys
			std::vector<double> x1, y1;
			x1.push_back(0.05f);
			y1.push_back(0.f);

			x1.push_back(0.1f);
			y1.push_back(0.1f);

			x1.push_back(0.325f);
			y1.push_back(0.2f);

			x1.push_back(0.5f);
			y1.push_back(0.21f);

			x1.push_back(0.65f);
			y1.push_back(0.55f);

			x1.push_back(0.9f);
			y1.push_back(0.61f);

			x1.push_back(1.f);
			y1.push_back(0.65f);
			tk::spline PVSpl;
			PVSpl.set_points(x1, y1);

			//Erosion
			std::vector<double> x2, y2;
			x2.push_back(0);
			y2.push_back(1);

			x2.push_back(0.15f);
			y2.push_back(0.55f);

			x2.push_back(0.3f);
			y2.push_back(0.41f);

			x2.push_back(0.32f);
			y2.push_back(0.45f);

			x2.push_back(0.45f);
			y2.push_back(PVSpl.solve(PV)[0]);

			x2.push_back(0.6f);
			y2.push_back(0.11f);

			x2.push_back(0.7f);
			y2.push_back(0.11f);

			x2.push_back(0.72f);
			y2.push_back(0.28f);

			x2.push_back(0.78f);
			y2.push_back(0.28f);

			x2.push_back(0.79f);
			y2.push_back(0.11f);

			x2.push_back(0.85f);
			y2.push_back(0.05f);

			x2.push_back(1.f);
			y2.push_back(0.05f);
			tk::spline ErosSpl;
			ErosSpl.set_points(x2, y2);

			//Continetalness
			std::vector<double> x3, y3;
			x3.push_back(0);
			y3.push_back(1);

			x3.push_back(0.09f);
			y3.push_back(0.025f);
			
			x3.push_back(0.345f);
			y3.push_back(0.025f);
			
			x3.push_back(0.385f);
			y3.push_back(0.435f);

			x3.push_back(0.515f);
			y3.push_back(0.435f);

			x3.push_back(0.525f);
			y3.push_back(0.9f);

			x3.push_back(0.55f);
			y3.push_back(0.915f);

			x3.push_back(0.72f);
			y3.push_back(ErosSpl.solve(Eros)[0]);

			x3.push_back(1.f);
			y3.push_back(0.98f);

			tk::spline ContSpl;
			ContSpl.set_points(x3, y3);
			
			worldHeight = ContSpl.solve(Cont)[0] * 100;

			return worldHeight;
		}

typo

Hi,

very nice code. Looks like you have a typo for case 2 below. Shouldn't it be 2.0*m_b0 instead of 2.0*m_b0*h?

if(x<m_x[0]) {
        // extrapolation to the left
        switch(order) {
        case 1:
            interpol=2.0*m_b0*h + m_c0;
            break;
        case 2:
            interpol=2.0*m_b0*h;
            break;
        default:
            interpol=0.0;
            break;
        }