+/**
+ * Get top most *point* on Bezier curve
+ */
+pair<BReal,BReal> Bezier::GetTop() const
+{
+ pair<BReal, BReal> tsols = BezierTurningPoints(y0,y1,y2,y3);
+ BReal tx0; BReal ty0;
+ BReal tx1; BReal ty1;
+ Evaluate(tx0, ty0, tsols.first);
+ Evaluate(tx1, ty1, tsols.second);
+ vector<const BReal*> v(4);
+ v[0] = &y0;
+ v[1] = &y3;
+ v[2] = &ty0;
+ v[3] = &ty1;
+ sort(v.begin(), v.end(), CompBRealByPtr);
+ pair<BReal,BReal> result;
+ result.second = *v[0];
+ if (v[0] == &y0)
+ {
+ result.first = x0;
+ }
+ else if (v[0] == &y3)
+ {
+ result.first = x3;
+ }
+ else if (v[0] == &ty0)
+ {
+ result.first = tx0;
+ }
+ else if (v[0] == &ty1)
+ {
+ result.first = tx1;
+ }
+ return result;
+}
+
+/**
+ * Get bottom most *point* on Bezier curve
+ */
+pair<BReal,BReal> Bezier::GetBottom() const
+{
+ pair<BReal, BReal> tsols = BezierTurningPoints(y0,y1,y2,y3);
+ BReal tx0; BReal ty0;
+ BReal tx1; BReal ty1;
+ Evaluate(tx0, ty0, tsols.first);
+ Evaluate(tx1, ty1, tsols.second);
+ vector<const BReal*> v(4);
+ v[0] = &y0;
+ v[1] = &y3;
+ v[2] = &ty0;
+ v[3] = &ty1;
+ sort(v.begin(), v.end(), CompBRealByPtr);
+ pair<BReal,BReal> result;
+ result.second = *v[3];
+ if (v[3] == &y0)
+ {
+ result.first = x0;
+ }
+ else if (v[3] == &y3)
+ {
+ result.first = x3;
+ }
+ else if (v[3] == &ty0)
+ {
+ result.first = tx0;
+ }
+ else if (v[3] == &ty1)
+ {
+ result.first = tx1;
+ }
+ return result;
+}
+
+/**
+ * Get left most *point* on Bezier curve
+ */
+pair<BReal,BReal> Bezier::GetLeft() const
+{
+ pair<BReal, BReal> tsols = BezierTurningPoints(x0,x1,x2,x3);
+ BReal tx0; BReal ty0;
+ BReal tx1; BReal ty1;
+ Evaluate(tx0, ty0, tsols.first);
+ Evaluate(tx1, ty1, tsols.second);
+ vector<const BReal*> v(4);
+ v[0] = &x0;
+ v[1] = &x3;
+ v[2] = &tx0;
+ v[3] = &tx1;
+ sort(v.begin(), v.end(), CompBRealByPtr);
+ pair<BReal,BReal> result;
+ result.first = *v[0];
+ if (v[0] == &x0)
+ {
+ result.second = y0;
+ }
+ else if (v[0] == &x3)
+ {
+ result.second = y3;
+ }
+ else if (v[0] == &tx0)
+ {
+ result.second = ty0;
+ }
+ else if (v[0] == &tx1)
+ {
+ result.second = ty1;
+ }
+ return result;
+}
+
+
+/**
+ * Get left most *point* on Bezier curve
+ */
+pair<BReal,BReal> Bezier::GetRight() const
+{
+ pair<BReal, BReal> tsols = BezierTurningPoints(x0,x1,x2,x3);
+ BReal tx0; BReal ty0;
+ BReal tx1; BReal ty1;
+ Evaluate(tx0, ty0, tsols.first);
+ Evaluate(tx1, ty1, tsols.second);
+ vector<const BReal*> v(4);
+ v[0] = &x0;
+ v[1] = &x3;
+ v[2] = &tx0;
+ v[3] = &tx1;
+ sort(v.begin(), v.end(), CompBRealByPtr);
+ pair<BReal,BReal> result;
+ result.first = *v[3];
+ if (v[3] == &x0)
+ {
+ result.second = y0;
+ }
+ else if (v[3] == &x3)
+ {
+ result.second = y3;
+ }
+ else if (v[3] == &tx0)
+ {
+ result.second = ty0;
+ }
+ else if (v[3] == &tx1)
+ {
+ result.second = ty1;
+ }
+ return result;
+}
+
+vector<BReal> Bezier::SolveXParam(const BReal & x) const
+{
+ BReal d(x0 - x);
+ BReal c((x1 - x0)*BReal(3));
+ BReal b((x2 - x1)*BReal(3) - c);
+ BReal a(x3 -x0 - c - b);
+ vector<BReal> results(SolveCubic(a, b, c, d));
+ for (unsigned i = 0; i < results.size(); ++i)
+ {
+ Vec2 p;
+ Evaluate(p.x, p.y, results[i]);
+ }
+ return results;
+}
+
+
+vector<BReal> Bezier::SolveYParam(const BReal & y) const
+{
+ BReal d(y0 - y);
+ BReal c((y1 - y0)*BReal(3));
+ BReal b((y2 - y1)*BReal(3) - c);
+ BReal a(y3 -y0 - c - b);
+ vector<BReal> results(SolveCubic(a, b, c, d));
+ for (unsigned i = 0; i < results.size(); ++i)
+ {
+ Vec2 p;
+ Evaluate(p.x, p.y, results[i]);
+ }
+ return results;
+}
+
+vector<Vec2> Bezier::Evaluate(const vector<BReal> & u) const
+{
+ vector<Vec2> result(u.size());
+ for (unsigned i = 0; i < u.size(); ++i)
+ {
+ Evaluate(result[i].x, result[i].y, u[i]);
+ }
+ return result;
+}
+