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);
103 for(int i = 0; i <= 100; ++i)
107 Real y = a*(x*x*x) + b*(x*x) + c*x + d;
108 if (((y < Real(0)) && (prevRes > Real(0))) || ((y > Real(0)) && (prevRes < Real(0))))
110 //Debug("Found root of %fx^3 + %fx^2 + %fx + %f at %f (%f)", a, b, c, d, x, y);
121 * Use dynamic programming / recursion
125 static unordered_map<int, int> dp;
126 static bool init = false;
132 auto it = dp.find(n);
135 int result = n*Factorial(n-1);
141 * Binomial coefficients
143 int BinomialCoeff(int n, int k)
145 return Factorial(n) / (Factorial(k) * Factorial(n-k));
149 * Bernstein Basis Polynomial
151 Real Bernstein(int k, int n, const Real & u)
153 return Real(BinomialCoeff(n, k)) * Power(u, k) * Power(Real(1.0) - u, n-k);
158 * Returns the parametric parameter at the turning point(s)
159 * In one coordinate direction
162 pair<Real, Real> BezierTurningPoints(const Real & p0, const Real & p1, const Real & p2, const Real & p3)
165 if (p1 == p2 && p2 == p3)
167 return pair<Real,Real>(0, 1);
169 Real a = (3*(p1-p2) + p3 - p0);
170 Real b = 2*(p2 - 2*p1 + p0);
175 return pair<Real, Real>(0,1);
179 return pair<Real, Real>(t, t);
181 //Debug("a, b, c are %f, %f, %f", Float(a), Float(b), Float(c));
184 //Debug("No real roots");
185 return pair<Real, Real>(0,1);
187 vector<Real> tsols = SolveQuadratic(a, b, c);
188 if (tsols.size() == 1)
189 return pair<Real,Real>(tsols[0], tsols[0]);
190 else if (tsols.size() == 0)
191 return pair<Real, Real>(0,1);
193 return pair<Real,Real>(tsols[0], tsols[1]);
197 inline bool CompRealByPtr(const Real * a, const Real * b)
203 * Get top most *point* on Bezier curve
205 pair<Real,Real> Bezier::GetTop() const
207 pair<Real, Real> tsols = BezierTurningPoints(y0,y1,y2,y3);
210 Evaluate(tx0, ty0, tsols.first);
211 Evaluate(tx1, ty1, tsols.second);
212 vector<const Real*> v(4);
217 sort(v.begin(), v.end(), CompRealByPtr);
218 pair<Real,Real> result;
219 result.second = *v[0];
224 else if (v[0] == &y3)
228 else if (v[0] == &ty0)
232 else if (v[0] == &ty1)
240 * Get bottom most *point* on Bezier curve
242 pair<Real,Real> Bezier::GetBottom() const
244 pair<Real, Real> tsols = BezierTurningPoints(y0,y1,y2,y3);
247 Evaluate(tx0, ty0, tsols.first);
248 Evaluate(tx1, ty1, tsols.second);
249 vector<const Real*> v(4);
254 sort(v.begin(), v.end(), CompRealByPtr);
255 pair<Real,Real> result;
256 result.second = *v[3];
261 else if (v[3] == &y3)
265 else if (v[3] == &ty0)
269 else if (v[3] == &ty1)
277 * Get left most *point* on Bezier curve
279 pair<Real,Real> Bezier::GetLeft() const
281 pair<Real, Real> tsols = BezierTurningPoints(x0,x1,x2,x3);
284 Evaluate(tx0, ty0, tsols.first);
285 Evaluate(tx1, ty1, tsols.second);
286 vector<const Real*> v(4);
291 sort(v.begin(), v.end(), CompRealByPtr);
292 pair<Real,Real> result;
293 result.first = *v[0];
298 else if (v[0] == &x3)
302 else if (v[0] == &tx0)
306 else if (v[0] == &tx1)
315 * Get left most *point* on Bezier curve
317 pair<Real,Real> Bezier::GetRight() const
319 pair<Real, Real> tsols = BezierTurningPoints(x0,x1,x2,x3);
322 Evaluate(tx0, ty0, tsols.first);
323 Evaluate(tx1, ty1, tsols.second);
324 vector<const Real*> v(4);
329 sort(v.begin(), v.end(), CompRealByPtr);
330 pair<Real,Real> result;
331 result.first = *v[3];
336 else if (v[3] == &x3)
340 else if (v[3] == &tx0)
344 else if (v[3] == &tx1)
351 vector<Real> Bezier::SolveXParam(const Real & x) const
354 Real c((x1 - x0)*Real(3));
355 Real b((x2 - x1)*Real(3) - c);
356 Real a(x3 -x0 - c - b);
357 vector<Real> results(SolveCubic(a, b, c, d));
358 for (unsigned i = 0; i < results.size(); ++i)
361 Evaluate(p.x, p.y, results[i]);
367 vector<Real> Bezier::SolveYParam(const Real & y) const
370 Real c((y1 - y0)*Real(3));
371 Real b((y2 - y1)*Real(3) - c);
372 Real a(y3 -y0 - c - b);
373 vector<Real> results(SolveCubic(a, b, c, d));
374 for (unsigned i = 0; i < results.size(); ++i)
377 Evaluate(p.x, p.y, results[i]);
382 vector<Vec2> Bezier::Evaluate(const vector<Real> & u) const
384 vector<Vec2> result(u.size());
385 for (unsigned i = 0; i < u.size(); ++i)
387 Evaluate(result[i].x, result[i].y, u[i]);
393 * Get Bounds Rectangle of Bezier
395 Rect Bezier::SolveBounds() const
398 pair<Real, Real> tsols = BezierTurningPoints(x0, x1, x2, x3);
400 Real tp0; Real tp1; Real o;
401 Evaluate(tp0, o, tsols.first);
402 Evaluate(tp1, o, tsols.second);
404 //Debug("x: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
406 vector<const Real*> v(4);
412 // Not using a lambda to keep this compiling on cabellera
413 sort(v.begin(), v.end(), CompRealByPtr);
416 result.w = *(v[3]) - result.x;
418 // Do the same thing for y component (wow this is a mess)
419 tsols = BezierTurningPoints(y0, y1, y2, y3);
420 Evaluate(o, tp0, tsols.first);
421 Evaluate(o, tp1, tsols.second);
424 //Debug("y: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
430 sort(v.begin(), v.end(), CompRealByPtr);
433 result.h = *(v[3]) - result.y;
435 //Debug("Solved Bezier %s bounds as %s", Str().c_str(), result.Str().c_str());