--- /dev/null
+#ifndef MSP_INTERPOLATE_BEZIERSPLINE_H_
+#define MSP_INTERPOLATE_BEZIERSPLINE_H_
+
+#include <vector>
+#include "spline.h"
+
+namespace Msp {
+namespace Interpolate {
+
+template<typename T, unsigned D>
+struct Bernstein
+{
+ static Polynomial<T, D> b(unsigned);
+};
+
+/**
+A spline composed of bezier curves. Each segment has D-1 control points
+between the knots which determine the shape of the curve. Bezier splines
+with 2- or 3-dimensional values are commonly used in graphics.
+*/
+template<typename T, unsigned D, unsigned N = 1>
+class BezierSpline: public Spline<T, D, N>
+{
+public:
+ using typename Spline<T, D, N>::Knot;
+
+ BezierSpline(const std::vector<Knot> &);
+};
+
+template<typename T, unsigned D>
+inline Polynomial<T, D> Bernstein<T, D>::b(unsigned k)
+{
+ if(k>D)
+ throw std::invalid_argument("Bernstein::b");
+
+ LinAl::Vector<T, D+1> coeff;
+
+ int c = ((D-k)%2 ? -1 : 1);
+ for(unsigned i=0; i<=D-k; ++i)
+ {
+ coeff[i] = c;
+ c = c*-1*static_cast<int>(D-k-i)/static_cast<int>(i+1);
+ }
+
+ unsigned n = 1;
+ for(unsigned i=0; i<k; ++i)
+ n *= D-i;
+ for(unsigned i=2; i<=k; ++i)
+ n /= i;
+
+ for(unsigned i=0; i<=D; ++i)
+ coeff[i] *= n;
+
+ return Polynomial<T, D>(coeff);
+}
+
+template<typename T, unsigned D, unsigned N>
+inline BezierSpline<T, D, N>::BezierSpline(const std::vector<Knot> &k):
+ Spline<T, D, N>(k.front())
+{
+ typedef SplineValue<T, N> SV;
+
+ if((k.size()-1)%D)
+ throw std::invalid_argument("BezierSpline::BezierSpline");
+
+ Polynomial<T, D> b[D+1];
+ for(unsigned i=0; i<=D; ++i)
+ b[i] = Bernstein<T, D>::b(i);
+
+ for(unsigned i=0; i+D<k.size(); i+=D)
+ {
+ T dx = k[i+D].x-k[i].x;
+ Polynomial<T, 1> t(LinAl::Vector<T, 2>(T(1)/dx, -k[i].x/dx));
+ Polynomial<T, D> p[N];
+ for(unsigned j=0; j<N; ++j)
+ {
+ for(unsigned c=0; c<=D; ++c)
+ p[j] += b[c]*SV::get(k[i+c].y, j);
+ p[j] = p[j](t);
+ }
+ this->add_segment(p, k[i+D].x);
+ }
+}
+
+} // namespace Interpolate
+} // namespace Msp
+
+#endif