bernstein polynoms for higher order bezier eval

include/spline/bernstein.h

+/**
+* \file bezier_curve.h
+* \brief class allowing to create a Bezier curve of dimension 1 <= n <= 3.
+* \author Steve T.
+* \version 0.1
+* \date 06/17/2013
+*/
+
+
+#ifndef _CLASS_BERNSTEIN
+#define _CLASS_BERNSTEIN
+
+#include "curve_abc.h"
+
+#include "MathDefs.h"
+
+#include <math.h>
+#include <vector>
+#include <stdexcept>
+
+namespace spline
+{
+///
+/// \brief Computes factorial of a number
+///
+unsigned int fact(const unsigned int n)
+{
+    assert(n>=0);
+    int res = 1;
+    for (int i=2 ; i <= n ; ++i)
+       res *= i;
+    return res;
+}
+
+///
+/// \brief Computes a binomal coefficient
+///
+unsigned int bin(const unsigned  int n, const unsigned  int k)
+{
+    return fact(n) / (fact(k) * fact(n - k));
+}
+
+/// \class Bernstein
+/// \brief Computes a Bernstein polynome
+///
+template <typename Numeric = double>
+struct Bern{
+Bern(const unsigned int m, const unsigned int i)
+    :m_minus_i(m - i)
+    ,i_(i)
+    ,bin_m_i_(bin(m,i)) {}
+
+~Bern(){}
+
+Numeric operator()(const Numeric u) const
+{
+    assert(u >= 0. && u <= 1.);
+    return bin_m_i_*(pow(u, i_)) *pow((1-u),m_minus_i);
+}
+
+Numeric m_minus_i;
+Numeric i_;
+Numeric bin_m_i_;
+};
+
+
+///
+/// \brief Computes all Bernstein polynomes for a certain degree
+///
+template <typename Numeric>
+std::vector<Bern<Numeric> > makeBernstein(const unsigned int n)
+{
+    std::vector<Bern<Numeric> > res;
+    for(unsigned int i = 0; i<= n; ++i)
+        res.push_back(Bern<Numeric>(n, i));
+    return res;
+}
+} // namespace spline
+#endif //_CLASS_BERNSTEIN
+

include/spline/bezier_curve.h

@@ -11,10 +11,12 @@
 #define _CLASS_BEZIERCURVE
 
 #include "curve_abc.h"
+#include "bernstein.h"
 
 #include "MathDefs.h"
 
 #include <vector>
+#include <stdexcept>
 
 namespace spline
 {
@@ -27,7 +29,8 @@ struct bezier_curve : public  curve_abc<Time, Numeric, Dim, Safe, Point>
 {
 	typedef Point 	point_t;
 	typedef Time 	time_t;
-	typedef Numeric	num_t;
+    typedef Numeric	num_t;
+    typedef std::vector<point_t,Eigen::aligned_allocator<point_t> > t_point_t;
 
 /* Constructors - destructors */
 	public:
@@ -39,11 +42,13 @@ struct bezier_curve : public  curve_abc<Time, Numeric, Dim, Safe, Point>
 	: minBound_(minBound)
 	, maxBound_(maxBound)
 	, size_(std::distance(PointsBegin, PointsEnd))
+    , bernstein_(spline::makeBernstein<num_t>(size_-1))
 	{
+        assert(bernstein_.size() == size_);
 		In it(PointsBegin);
 		if(Safe && (size_<=1 || minBound == maxBound))
 		{
-            throw std::out_of_range("TODO"); // TODO
+            throw std::out_of_range("can't create bezier min bound is higher than max bound"); // TODO
 		}
 		for(; it != PointsEnd; ++it)
 		{
@@ -58,8 +63,8 @@ struct bezier_curve : public  curve_abc<Time, Numeric, Dim, Safe, Point>
 	}
 
 	private:
-	bezier_curve(const bezier_curve&);
-	bezier_curve& operator=(const bezier_curve&);
+//	bezier_curve(const bezier_curve&);
+//  bezier_curve& operator=(const bezier_curve&);
 /* Constructors - destructors */
 
 /*Operations*/
@@ -72,7 +77,7 @@ struct bezier_curve : public  curve_abc<Time, Numeric, Dim, Safe, Point>
 		num_t nT = (t - minBound_) / (maxBound_ - minBound_);
 		if(Safe &! (0 <= nT && nT <= 1))
 		{
-            throw std::out_of_range("TODO"); // TODO
+            throw std::out_of_range("can't evaluate bezier curve, out of range"); // TODO
         }
 		else
 		{
@@ -87,15 +92,33 @@ struct bezier_curve : public  curve_abc<Time, Numeric, Dim, Safe, Point>
                        				+ 2 * pts_[1] * nT * dt
 						+ pts_[2] * nT * nT;
 				break;
-				default :
+                case 4 :
 					return 	pts_[0] * dt * dt * dt
 						+ 3 * pts_[1] * nT * dt * dt 
 						+ 3 * pts_[2] * nT * nT * dt 
 						+ pts_[3] * nT * nT *nT;
+                default :
+                    return evalBernstein(dt);
 				break;
 			}
 		}
 	}
+
+    ///
+    /// \brief Evaluates all Bernstein polynomes for a certain degree
+    ///
+    point_t evalBernstein(const Numeric u) const
+    {
+        point_t res = point_t::Zero();
+        typename t_point_t::const_iterator pts_it = pts_.begin();
+        for(typename std::vector<Bern<Numeric> >::const_iterator cit = bernstein_.begin();
+            cit !=bernstein_.end(); ++cit, ++pts_it)
+        {
+            res += cit->operator()(u) * (*pts_it);
+        }
+        return res;
+    }
+
 /*Operations*/
 
 /*Helpers*/
@@ -104,12 +127,14 @@ struct bezier_curve : public  curve_abc<Time, Numeric, Dim, Safe, Point>
 /*Helpers*/
 
 	public:
-	const int size_;
-	const time_t minBound_, maxBound_;
+    const time_t minBound_, maxBound_;
+    const std::size_t size_;
+    const std::vector<Bern<Numeric> > bernstein_;
 	
-	private:
-	typedef std::vector<Point,Eigen::aligned_allocator<Point> > T_Vector;
-	T_Vector  pts_;
+    private:
+    t_point_t  pts_;
+
+    //storing bernstein polynoms, even in low dimension
 };
 }
 #endif //_CLASS_BEZIERCURVE

src/tests/spline_test/Main.cpp

@@ -181,30 +181,40 @@ void BezierCurveTest(bool& error)
 	res1 = cf4(2); 
     ComparePoints(d, res1, errMsg + "3(1) ", error);
 
+    //testing bernstein polynomes
+    std::string errMsg2("In test BezierCurveTest ; Bernstein polynoms do not evaluate as analytical evaluation");
+    for(double d = 0.; d <1.; d+=0.1)
+    {
+        ComparePoints( cf3.evalBernstein(d) , cf3 (d), errMsg2, error);
+    }
+
+    bool error_in(true);
 	try
 	{
 		cf(-0.4);
 	}
 	catch(...)
 	{
-		error = false;
+        error_in = false;
 	}
-	if(error)
+    if(error_in)
 	{
 		std::cout << "Evaluation of bezier cf error, -0.4 should be an out of range value\n";
+        error = true;
 	}
-	error = true;	
+    error_in = true;
 	try
 	{
 		cf(1.1);
 	}
 	catch(...)
 	{
-		error = false;
+        error_in = false;
 	}
-	if(error)
+    if(error_in)
 	{
 		std::cout << "Evaluation of bezier cf error, 1.1 should be an out of range value\n";
+        error = true;
 	}
 	if(cf.max() != 1)
 	{