Add a bezier spline type

[libs/math.git] / source / interpolate / bezierspline.h

diff --git a/source/interpolate/bezierspline.h b/source/interpolate/bezierspline.h
new file mode 100644 (file)
index 0000000..c6e295a
--- /dev/null
+++ b/source/interpolate/bezierspline.h
@@ -0,0 +1,88 @@
+#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