3 #include <unordered_map>
12 vector<Real> SolveQuadratic(const Real & a, const Real & b, const Real & c, const Real & min, const Real & max)
14 vector<Real> roots; roots.reserve(2);
17 roots.push_back(-c/b);
20 Real disc(b*b - Real(4)*a*c);
28 if (x >= min && x <= max)
33 Real x0((-b - Sqrt(b*b - Real(4)*a*c))/(Real(2)*a));
34 Real x1((-b + Sqrt(b*b - Real(4)*a*c))/(Real(2)*a));
41 if (x0 >= min && x0 <= max)
43 if (x1 >= min && x1 <= max)
49 * Finds the root (if it exists) in a monotonicly in(de)creasing segment of a Cubic
52 static void CubicSolveSegment(vector<Real> & roots, const Real & a, const Real & b, const Real & c, const Real & d, Real & tl, Real & tu, const Real & delta)
54 Real l = a*tl*tl*tl + b*tl*tl + c*tl + d;
55 Real u = a*tu*tu*tu + b*tu*tu + c*tu + d;
56 if ((l < 0 && u < 0) || (l > 0 && u > 0))
59 bool negative = (u < l); // lower point > 0, upper point < 0
60 while (tu - tl > delta)
64 Real m = a*t*t*t + b*t*t + c*t + d;
80 //Debug("Delta is %f (%f - %f -> %f)", tu-tl, tu, tl, t);
84 vector<Real> SolveCubic(const Real & a, const Real & b, const Real & c, const Real & d, const Real & min, const Real & max, const Real & delta)
86 vector<Real> roots; roots.reserve(3);
89 vector<Real> turns(SolveQuadratic(a*3, b*2, c));
90 //Debug("%u turning points", turns.size());
91 for (unsigned i = 1; i < turns.size(); ++i)
94 CubicSolveSegment(roots, a, b, c, d, tl, tu,delta);
98 CubicSolveSegment(roots, a, b, c, d, tl, tu,delta);
104 * Use dynamic programming / recursion
108 static unordered_map<int, int> dp;
109 static bool init = false;
115 auto it = dp.find(n);
118 int result = n*Factorial(n-1);
124 * Binomial coefficients
126 int BinomialCoeff(int n, int k)
128 return Factorial(n) / (Factorial(k) * Factorial(n-k));
132 * Bernstein Basis Polynomial
134 Real Bernstein(int k, int n, const Real & u)
136 return Real(BinomialCoeff(n, k)) * Power(u, k) * Power(Real(1.0) - u, n-k);
141 * Returns the parametric parameter at the turning point(s)
142 * In one coordinate direction
145 pair<Real, Real> BezierTurningPoints(const Real & p0, const Real & p1, const Real & p2, const Real & p3)
148 if (p1 == p2 && p2 == p3)
150 return pair<Real,Real>(0, 1);
152 Real a = ((p1-p2)*3 + p3 - p0);
153 Real b = (p2 - p1*2 + p0)*2;
158 return pair<Real, Real>(0,1);
162 return pair<Real, Real>(t, t);
164 //Debug("a, b, c are %f, %f, %f", Float(a), Float(b), Float(c));
167 //Debug("No real roots");
168 return pair<Real, Real>(0,1);
170 vector<Real> tsols = SolveQuadratic(a, b, c);
171 if (tsols.size() == 1)
172 return pair<Real,Real>(tsols[0], tsols[0]);
173 else if (tsols.size() == 0)
174 return pair<Real, Real>(0,1);
176 return pair<Real,Real>(tsols[0], tsols[1]);
180 inline bool CompRealByPtr(const Real * a, const Real * b)
186 * Get top most *point* on Bezier curve
188 pair<Real,Real> Bezier::GetTop() const
190 pair<Real, Real> tsols = BezierTurningPoints(y0,y1,y2,y3);
193 Evaluate(tx0, ty0, tsols.first);
194 Evaluate(tx1, ty1, tsols.second);
195 vector<const Real*> v(4);
200 sort(v.begin(), v.end(), CompRealByPtr);
201 pair<Real,Real> result;
202 result.second = *v[0];
207 else if (v[0] == &y3)
211 else if (v[0] == &ty0)
215 else if (v[0] == &ty1)
223 * Get bottom most *point* on Bezier curve
225 pair<Real,Real> Bezier::GetBottom() const
227 pair<Real, Real> tsols = BezierTurningPoints(y0,y1,y2,y3);
230 Evaluate(tx0, ty0, tsols.first);
231 Evaluate(tx1, ty1, tsols.second);
232 vector<const Real*> v(4);
237 sort(v.begin(), v.end(), CompRealByPtr);
238 pair<Real,Real> result;
239 result.second = *v[3];
244 else if (v[3] == &y3)
248 else if (v[3] == &ty0)
252 else if (v[3] == &ty1)
260 * Get left most *point* on Bezier curve
262 pair<Real,Real> Bezier::GetLeft() const
264 pair<Real, Real> tsols = BezierTurningPoints(x0,x1,x2,x3);
267 Evaluate(tx0, ty0, tsols.first);
268 Evaluate(tx1, ty1, tsols.second);
269 vector<const Real*> v(4);
274 sort(v.begin(), v.end(), CompRealByPtr);
275 pair<Real,Real> result;
276 result.first = *v[0];
281 else if (v[0] == &x3)
285 else if (v[0] == &tx0)
289 else if (v[0] == &tx1)
298 * Get left most *point* on Bezier curve
300 pair<Real,Real> Bezier::GetRight() const
302 pair<Real, Real> tsols = BezierTurningPoints(x0,x1,x2,x3);
305 Evaluate(tx0, ty0, tsols.first);
306 Evaluate(tx1, ty1, tsols.second);
307 vector<const Real*> v(4);
312 sort(v.begin(), v.end(), CompRealByPtr);
313 pair<Real,Real> result;
314 result.first = *v[3];
319 else if (v[3] == &x3)
323 else if (v[3] == &tx0)
327 else if (v[3] == &tx1)
334 vector<Real> Bezier::SolveXParam(const Real & x) const
337 Real c((x1 - x0)*Real(3));
338 Real b((x2 - x1)*Real(3) - c);
339 Real a(x3 -x0 - c - b);
340 vector<Real> results(SolveCubic(a, b, c, d));
341 for (unsigned i = 0; i < results.size(); ++i)
344 Evaluate(p.x, p.y, results[i]);
350 vector<Real> Bezier::SolveYParam(const Real & y) const
353 Real c((y1 - y0)*Real(3));
354 Real b((y2 - y1)*Real(3) - c);
355 Real a(y3 -y0 - c - b);
356 vector<Real> results(SolveCubic(a, b, c, d));
357 for (unsigned i = 0; i < results.size(); ++i)
360 Evaluate(p.x, p.y, results[i]);
365 vector<Vec2> Bezier::Evaluate(const vector<Real> & u) const
367 vector<Vec2> result(u.size());
368 for (unsigned i = 0; i < u.size(); ++i)
370 Evaluate(result[i].x, result[i].y, u[i]);
376 * Get Bounds Rectangle of Bezier
378 Rect Bezier::SolveBounds() const
381 pair<Real, Real> tsols = BezierTurningPoints(x0, x1, x2, x3);
383 Real tp0; Real tp1; Real o;
384 Evaluate(tp0, o, tsols.first);
385 Evaluate(tp1, o, tsols.second);
387 //Debug("x: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
389 vector<const Real*> v(4);
395 // Not using a lambda to keep this compiling on cabellera
396 sort(v.begin(), v.end(), CompRealByPtr);
399 result.w = *(v[3]) - result.x;
401 // Do the same thing for y component (wow this is a mess)
402 tsols = BezierTurningPoints(y0, y1, y2, y3);
403 Evaluate(o, tp0, tsols.first);
404 Evaluate(o, tp1, tsols.second);
407 //Debug("y: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
413 sort(v.begin(), v.end(), CompRealByPtr);
416 result.h = *(v[3]) - result.y;
418 //Debug("Solved Bezier %s bounds as %s", Str().c_str(), result.Str().c_str());