Cubic spline interpolation is a simple way to obtain a smooth curve from a set of discrete points (knots). It has both C1 (first derivative) and C2 (second derivative) continuity, enabling it to produce a continuous piecewise function given a set of data points.

From the algorithm detailed here I have implemented a clamped cubic spline class in C++. It is templated on the type of X and Y, allowing for use of scalar or vector types. It requires only that X and Y define the operators +, -, *, and /, and that Y have a constructor that takes a single scalar argument, initializing all elements to the supplied value.

Usage is very simple. For example, to interpolate a 2D location over time, try this:

#import "glm.hpp" #import <vector> std::vector<float> times; std::vector<glm::vec2> points; times.push_back(0); points.push_back(glm::vec2(0)); times.push_back(.5f); points.push_back(glm::vec2(-303, -572)); times.push_back(1); points.push_back(glm::vec2(-250, -980)); /* Create the spline interpolating the position over time */ Spline<float, glm::vec2> sp(times, points); /* Compute the interpolated position at time 0.75f */ glm::vec2 value(sp[.75f]);

The code for the spline class is below:

/* "THE BEER-WARE LICENSE" (Revision 42): Devin Lane wrote this file. As long as you retain * this notice you can do whatever you want with this stuff. If we meet some day, and you * think this stuff is worth it, you can buy me a beer in return. */ #include <vector> #include <iostream> /** Templated on type of X, Y. X and Y must have operator +, -, *, /. Y must have defined * a constructor that takes a scalar. */ template <typename X, typename Y> class Spline { public: /** An empty, invalid spline */ Spline() {} /** A spline with x and y values */ Spline(const std::vector<X>& x, const std::vector<Y>& y) { if (x.size() != y.size()) { std::cerr << "X and Y must be the same size " << std::endl; return; } if (x.size() < 3) { std::cerr << "Must have at least three points for interpolation" << std::endl; return; } typedef typename std::vector<X>::difference_type size_type; size_type n = y.size() - 1; std::vector<Y> b(n), d(n), a(n), c(n+1), l(n+1), u(n+1), z(n+1); std::vector<X> h(n+1); l[0] = Y(1); u[0] = Y(0); z[0] = Y(0); h[0] = x[1] - x[0]; for (size_type i = 1; i < n; i++) { h[i] = x[i+1] - x[i]; l[i] = Y(2 * (x[i+1] - x[i-1])) - Y(h[i-1]) * u[i-1]; u[i] = Y(h[i]) / l[i]; a[i] = (Y(3) / Y(h[i])) * (y[i+1] - y[i]) - (Y(3) / Y(h[i-1])) * (y[i] - y[i-1]); z[i] = (a[i] - Y(h[i-1]) * z[i-1]) / l[i]; } l[n] = Y(1); z[n] = c[n] = Y(0); for (size_type j = n-1; j >= 0; j--) { c[j] = z[j] - u[j] * c[j+1]; b[j] = (y[j+1] - y[j]) / Y(h[j]) - (Y(h[j]) * (c[j+1] + Y(2) * c[j])) / Y(3); d[j] = (c[j+1] - c[j]) / Y(3 * h[j]); } for (size_type i = 0; i < n; i++) { mElements.push_back(Element(x[i], y[i], b[i], c[i], d[i])); } } virtual ~Spline() {} Y operator[](const X& x) const { return interpolate(x); } Y interpolate(const X&x) const { if (mElements.size() == 0) return Y(); typename std::vector<element_type>::const_iterator it; it = std::lower_bound(mElements.begin(), mElements.end(), element_type(x)); if (it != mElements.begin()) { it--; } return it->eval(x); } std::vector<Y> operator[](const std::vector<X>& xx) const { return interpolate(xx); } /* Evaluate at multiple locations, assuming xx is sorted ascending */ std::vector<Y> interpolate(const std::vector<X>& xx) const { if (mElements.size() == 0) return std::vector<Y>(xx.size()); typename std::vector<X>::const_iterator it; typename std::vector<element_type>::const_iterator it2; it2 = mElements.begin(); std::vector<Y> ys; for (it = xx.begin(); it != xx.end(); it++) { it2 = std::lower_bound(it2, mElements.end(), element_type(*it)); if (it2 != mElements.begin()) { it2--; } ys.push_back(it2->eval(*it)); } return ys; } protected: class Element { public: Element(X _x) : x(_x) {} Element(X _x, Y _a, Y _b, Y _c, Y _d) : x(_x), a(_a), b(_b), c(_c), d(_d) {} Y eval(const X& xx) const { X xix(xx - x); return a + b * xix + c * (xix * xix) + d * (xix * xix * xix); } bool operator<(const Element& e) const { return x < e.x; } bool operator<(const X& xx) const { return x < xx; } X x; Y a, b, c, d; }; typedef Element element_type; std::vector<element_type> mElements; };

How would I get results out of this to use them to draw a graph on a UIView in Objective-C?

Thank you so much for posting such code.

This will be so much helpful.

Thanks for posting this!

Finally after more than one month I’ve finally been able to implement the CubicSpline interpolation!!!

Thank you for the code !

It works great

How one can use that class. Could any body help to provide the main code. regards

Thanks. Is there any usage example?

Thanks Man! This is the best spline library I saw today and I saw quite a few. I’ll definitely buy you a beer!

Thanks. Good stuff. I’m trying out a modified version of your library that add support for calculating bezier control points. You can find it at github (Spline directory of Engauge6). Some parts may be useful to you, although I’m not recommending you remove the templates like I did for my very specific use case. BTW, any suggestions for delivering on the beer fee?