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))
57 Debug("Discarding segment (no roots) l = %f (%f), u = %f (%f)", Double(tl), Double(l), Double(tu), Double(u));
60 bool negative = (u < l); // lower point > 0, upper point < 0
61 Debug("%ft^3 + %ft^2 + %ft + %f is negative (%f < %f) %d", Double(a),Double(b),Double(c),Double(d),Double(u),Double(l), negative);
62 while (tu - tl > delta)
66 Real m = a*t*t*t + b*t*t + c*t + d;
82 //Debug("Delta is %f (%f - %f -> %f)", tu-tl, tu, tl, t);
86 vector<Real> SolveCubic(const Real & a, const Real & b, const Real & c, const Real & d, const Real & min, const Real & max, const Real & delta)
88 vector<Real> roots; roots.reserve(3);
91 vector<Real> turns(SolveQuadratic(a*3, b*2, c));
92 Debug("%u turning points", turns.size());
93 for (unsigned i = 1; i < turns.size(); ++i)
96 CubicSolveSegment(roots, a, b, c, d, tl, tu,delta);
100 CubicSolveSegment(roots, a, b, c, d, tl, tu,delta);
106 * Use dynamic programming / recursion
110 static unordered_map<int, int> dp;
111 static bool init = false;
117 auto it = dp.find(n);
120 int result = n*Factorial(n-1);
126 * Binomial coefficients
128 int BinomialCoeff(int n, int k)
130 return Factorial(n) / (Factorial(k) * Factorial(n-k));
134 * Bernstein Basis Polynomial
136 Real Bernstein(int k, int n, const Real & u)
138 return Real(BinomialCoeff(n, k)) * Power(u, k) * Power(Real(1.0) - u, n-k);
143 * Returns the parametric parameter at the turning point(s)
144 * In one coordinate direction
147 pair<Real, Real> BezierTurningPoints(const Real & p0, const Real & p1, const Real & p2, const Real & p3)
150 if (p1 == p2 && p2 == p3)
152 return pair<Real,Real>(0, 1);
154 Real a = ((p1-p2)*3 + p3 - p0);
155 Real b = (p2 - p1*2 + p0)*2;
160 return pair<Real, Real>(0,1);
164 return pair<Real, Real>(t, t);
166 //Debug("a, b, c are %f, %f, %f", Float(a), Float(b), Float(c));
169 //Debug("No real roots");
170 return pair<Real, Real>(0,1);
172 vector<Real> tsols = SolveQuadratic(a, b, c);
173 if (tsols.size() == 1)
174 return pair<Real,Real>(tsols[0], tsols[0]);
175 else if (tsols.size() == 0)
176 return pair<Real, Real>(0,1);
178 return pair<Real,Real>(tsols[0], tsols[1]);
182 inline bool CompRealByPtr(const Real * a, const Real * b)
188 * Get top most *point* on Bezier curve
190 pair<Real,Real> Bezier::GetTop() const
192 pair<Real, Real> tsols = BezierTurningPoints(y0,y1,y2,y3);
195 Evaluate(tx0, ty0, tsols.first);
196 Evaluate(tx1, ty1, tsols.second);
197 vector<const Real*> v(4);
202 sort(v.begin(), v.end(), CompRealByPtr);
203 pair<Real,Real> result;
204 result.second = *v[0];
209 else if (v[0] == &y3)
213 else if (v[0] == &ty0)
217 else if (v[0] == &ty1)
225 * Get bottom most *point* on Bezier curve
227 pair<Real,Real> Bezier::GetBottom() const
229 pair<Real, Real> tsols = BezierTurningPoints(y0,y1,y2,y3);
232 Evaluate(tx0, ty0, tsols.first);
233 Evaluate(tx1, ty1, tsols.second);
234 vector<const Real*> v(4);
239 sort(v.begin(), v.end(), CompRealByPtr);
240 pair<Real,Real> result;
241 result.second = *v[3];
246 else if (v[3] == &y3)
250 else if (v[3] == &ty0)
254 else if (v[3] == &ty1)
262 * Get left most *point* on Bezier curve
264 pair<Real,Real> Bezier::GetLeft() const
266 pair<Real, Real> tsols = BezierTurningPoints(x0,x1,x2,x3);
269 Evaluate(tx0, ty0, tsols.first);
270 Evaluate(tx1, ty1, tsols.second);
271 vector<const Real*> v(4);
276 sort(v.begin(), v.end(), CompRealByPtr);
277 pair<Real,Real> result;
278 result.first = *v[0];
283 else if (v[0] == &x3)
287 else if (v[0] == &tx0)
291 else if (v[0] == &tx1)
300 * Get left most *point* on Bezier curve
302 pair<Real,Real> Bezier::GetRight() const
304 pair<Real, Real> tsols = BezierTurningPoints(x0,x1,x2,x3);
307 Evaluate(tx0, ty0, tsols.first);
308 Evaluate(tx1, ty1, tsols.second);
309 vector<const Real*> v(4);
314 sort(v.begin(), v.end(), CompRealByPtr);
315 pair<Real,Real> result;
316 result.first = *v[3];
321 else if (v[3] == &x3)
325 else if (v[3] == &tx0)
329 else if (v[3] == &tx1)
336 vector<Real> Bezier::SolveXParam(const Real & x) const
339 Real c((x1 - x0)*Real(3));
340 Real b((x2 - x1)*Real(3) - c);
341 Real a(x3 -x0 - c - b);
342 vector<Real> results(SolveCubic(a, b, c, d));
343 for (unsigned i = 0; i < results.size(); ++i)
346 Evaluate(p.x, p.y, results[i]);
352 vector<Real> Bezier::SolveYParam(const Real & y) const
355 Real c((y1 - y0)*Real(3));
356 Real b((y2 - y1)*Real(3) - c);
357 Real a(y3 -y0 - c - b);
358 vector<Real> results(SolveCubic(a, b, c, d));
359 for (unsigned i = 0; i < results.size(); ++i)
362 Evaluate(p.x, p.y, results[i]);
367 vector<Vec2> Bezier::Evaluate(const vector<Real> & u) const
369 vector<Vec2> result(u.size());
370 for (unsigned i = 0; i < u.size(); ++i)
372 Evaluate(result[i].x, result[i].y, u[i]);
378 * Get Bounds Rectangle of Bezier
380 Rect Bezier::SolveBounds() const
383 pair<Real, Real> tsols = BezierTurningPoints(x0, x1, x2, x3);
385 Real tp0; Real tp1; Real o;
386 Evaluate(tp0, o, tsols.first);
387 Evaluate(tp1, o, tsols.second);
389 //Debug("x: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
391 vector<const Real*> v(4);
397 // Not using a lambda to keep this compiling on cabellera
398 sort(v.begin(), v.end(), CompRealByPtr);
401 result.w = *(v[3]) - result.x;
403 // Do the same thing for y component (wow this is a mess)
404 tsols = BezierTurningPoints(y0, y1, y2, y3);
405 Evaluate(o, tp0, tsols.first);
406 Evaluate(o, tp1, tsols.second);
409 //Debug("y: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
415 sort(v.begin(), v.end(), CompRealByPtr);
418 result.h = *(v[3]) - result.y;
420 //Debug("Solved Bezier %s bounds as %s", Str().c_str(), result.Str().c_str());