3 #include <unordered_map>
14 vector<Real> SolveQuadratic(const Real & a, const Real & b, const Real & c, const Real & min, const Real & max)
16 vector<Real> roots; roots.reserve(2);
19 roots.push_back(-c/b);
22 Real disc(b*b - Real(4)*a*c);
30 if (x >= min && x <= max)
35 Real x0((-b - Sqrt(b*b - Real(4)*a*c))/(Real(2)*a));
36 Real x1((-b + Sqrt(b*b - Real(4)*a*c))/(Real(2)*a));
43 if (x0 >= min && x0 <= max)
45 if (x1 >= min && x1 <= max)
51 * Finds the root (if it exists) in a monotonicly in(de)creasing segment of a Cubic
54 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)
56 Real l = a*tl*tl*tl + b*tl*tl + c*tl + d;
57 Real u = a*tu*tu*tu + b*tu*tu + c*tu + d;
58 if ((l < 0 && u < 0) || (l > 0 && u > 0))
60 //Debug("Discarding segment (no roots) l = %f (%f), u = %f (%f)", Double(tl), Double(l), Double(tu), Double(u));
64 bool negative = (u < l); // lower point > 0, upper point < 0
65 //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);
66 while (tu - tl > delta)
70 Real m = a*t*t*t + b*t*t + c*t + d;
86 //Debug("Delta is %f (%f - %f -> %f)", tu-tl, tu, tl, t);
90 vector<Real> SolveCubic(const Real & a, const Real & b, const Real & c, const Real & d, const Real & min, const Real & max, const Real & delta)
92 vector<Real> roots; roots.reserve(3);
95 vector<Real> turns(SolveQuadratic(a*3, b*2, c));
96 //Debug("%u turning points", turns.size());
97 for (unsigned i = 1; i < turns.size(); ++i)
100 CubicSolveSegment(roots, a, b, c, d, tl, tu,delta);
104 CubicSolveSegment(roots, a, b, c, d, tl, tu,delta);
110 * Use dynamic programming / recursion
114 static unordered_map<int, int> dp;
115 static bool init = false;
121 auto it = dp.find(n);
124 int result = n*Factorial(n-1);
130 * Binomial coefficients
132 int BinomialCoeff(int n, int k)
134 return Factorial(n) / (Factorial(k) * Factorial(n-k));
138 * Bernstein Basis Polynomial
140 Real Bernstein(int k, int n, const Real & u)
142 return Real(BinomialCoeff(n, k)) * Power(u, k) * Power(Real(1.0) - u, n-k);
147 * Returns the parametric parameter at the turning point(s)
148 * In one coordinate direction
151 pair<Real, Real> BezierTurningPoints(const Real & p0, const Real & p1, const Real & p2, const Real & p3)
154 if (p1 == p2 && p2 == p3)
156 return pair<Real,Real>(0, 1);
158 Real a = ((p1-p2)*3 + p3 - p0);
159 Real b = (p2 - p1*2 + p0)*2;
164 return pair<Real, Real>(0,1);
168 return pair<Real, Real>(t, t);
170 //Debug("a, b, c are %f, %f, %f", Float(a), Float(b), Float(c));
173 //Debug("No real roots");
174 return pair<Real, Real>(0,1);
176 vector<Real> tsols = SolveQuadratic(a, b, c);
177 if (tsols.size() == 1)
178 return pair<Real,Real>(tsols[0], tsols[0]);
179 else if (tsols.size() == 0)
180 return pair<Real, Real>(0,1);
182 return pair<Real,Real>(tsols[0], tsols[1]);
186 inline bool CompRealByPtr(const Real * a, const Real * b)
192 * Get top most *point* on Bezier curve
194 pair<Real,Real> Bezier::GetTop() const
196 pair<Real, Real> tsols = BezierTurningPoints(y0,y1,y2,y3);
199 Evaluate(tx0, ty0, tsols.first);
200 Evaluate(tx1, ty1, tsols.second);
201 vector<const Real*> v(4);
206 sort(v.begin(), v.end(), CompRealByPtr);
207 pair<Real,Real> result;
208 result.second = *v[0];
213 else if (v[0] == &y3)
217 else if (v[0] == &ty0)
221 else if (v[0] == &ty1)
229 * Get bottom most *point* on Bezier curve
231 pair<Real,Real> Bezier::GetBottom() const
233 pair<Real, Real> tsols = BezierTurningPoints(y0,y1,y2,y3);
236 Evaluate(tx0, ty0, tsols.first);
237 Evaluate(tx1, ty1, tsols.second);
238 vector<const Real*> v(4);
243 sort(v.begin(), v.end(), CompRealByPtr);
244 pair<Real,Real> result;
245 result.second = *v[3];
250 else if (v[3] == &y3)
254 else if (v[3] == &ty0)
258 else if (v[3] == &ty1)
266 * Get left most *point* on Bezier curve
268 pair<Real,Real> Bezier::GetLeft() const
270 pair<Real, Real> tsols = BezierTurningPoints(x0,x1,x2,x3);
273 Evaluate(tx0, ty0, tsols.first);
274 Evaluate(tx1, ty1, tsols.second);
275 vector<const Real*> v(4);
280 sort(v.begin(), v.end(), CompRealByPtr);
281 pair<Real,Real> result;
282 result.first = *v[0];
287 else if (v[0] == &x3)
291 else if (v[0] == &tx0)
295 else if (v[0] == &tx1)
304 * Get left most *point* on Bezier curve
306 pair<Real,Real> Bezier::GetRight() const
308 pair<Real, Real> tsols = BezierTurningPoints(x0,x1,x2,x3);
311 Evaluate(tx0, ty0, tsols.first);
312 Evaluate(tx1, ty1, tsols.second);
313 vector<const Real*> v(4);
318 sort(v.begin(), v.end(), CompRealByPtr);
319 pair<Real,Real> result;
320 result.first = *v[3];
325 else if (v[3] == &x3)
329 else if (v[3] == &tx0)
333 else if (v[3] == &tx1)
340 vector<Real> Bezier::SolveXParam(const Real & x) const
343 Real c((x1 - x0)*Real(3));
344 Real b((x2 - x1)*Real(3) - c);
345 Real a(x3 -x0 - c - b);
346 vector<Real> results(SolveCubic(a, b, c, d));
347 for (unsigned i = 0; i < results.size(); ++i)
350 Evaluate(p.x, p.y, results[i]);
356 vector<Real> Bezier::SolveYParam(const Real & y) const
359 Real c((y1 - y0)*Real(3));
360 Real b((y2 - y1)*Real(3) - c);
361 Real a(y3 -y0 - c - b);
362 vector<Real> results(SolveCubic(a, b, c, d));
363 for (unsigned i = 0; i < results.size(); ++i)
366 Evaluate(p.x, p.y, results[i]);
371 vector<Vec2> Bezier::Evaluate(const vector<Real> & u) const
373 vector<Vec2> result(u.size());
374 for (unsigned i = 0; i < u.size(); ++i)
376 Evaluate(result[i].x, result[i].y, u[i]);
382 * Get Bounds Rectangle of Bezier
384 Rect Bezier::SolveBounds() const
387 pair<Real, Real> tsols = BezierTurningPoints(x0, x1, x2, x3);
389 Real tp0; Real tp1; Real o;
390 Evaluate(tp0, o, tsols.first);
391 Evaluate(tp1, o, tsols.second);
393 //Debug("x: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
395 vector<const Real*> v(4);
401 // Not using a lambda to keep this compiling on cabellera
402 sort(v.begin(), v.end(), CompRealByPtr);
405 result.w = *(v[3]) - result.x;
407 // Do the same thing for y component (wow this is a mess)
408 tsols = BezierTurningPoints(y0, y1, y2, y3);
409 Evaluate(o, tp0, tsols.first);
410 Evaluate(o, tp1, tsols.second);
413 //Debug("y: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
419 sort(v.begin(), v.end(), CompRealByPtr);
422 result.h = *(v[3]) - result.y;
424 //Debug("Solved Bezier %s bounds as %s", Str().c_str(), result.Str().c_str());