X-Git-Url: http://git.tdb.fi/?p=libs%2Fmath.git;a=blobdiff_plain;f=source%2Finterpolate%2Fbezierspline.h;fp=source%2Finterpolate%2Fbezierspline.h;h=c6e295a3dc40072aa9804de0977d280d181bfa19;hp=0000000000000000000000000000000000000000;hb=c854dd42c57b0430558baed5151449c83e24507f;hpb=92480ec2f41502e3cbdfe1fe3ca703c4521ff747

diff --git a/source/interpolate/bezierspline.h b/source/interpolate/bezierspline.h
new file mode 100644
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