3 #include <unordered_map>
14 vector<BReal> SolveQuadratic(const BReal & a, const BReal & b, const BReal & c, const BReal & min, const BReal & max)
16 vector<BReal> roots; roots.reserve(2);
19 roots.push_back(-c/b);
22 BReal disc(b*b - BReal(4)*a*c);
29 BReal x(-b/BReal(2)*a);
30 if (x >= min && x <= max)
35 BReal x0((-b - Sqrt(b*b - BReal(4)*a*c))/(BReal(2)*a));
36 BReal x1((-b + Sqrt(b*b - BReal(4)*a*c))/(BReal(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<BReal> & roots, const BReal & a, const BReal & b, const BReal & c, const BReal & d, BReal & tl, BReal & tu, const BReal & delta)
56 BReal l = a*tl*tl*tl + b*tl*tl + c*tl + d;
57 BReal 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 BReal 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<BReal> SolveCubic(const BReal & a, const BReal & b, const BReal & c, const BReal & d, const BReal & min, const BReal & max, const BReal & delta)
92 vector<BReal> roots; roots.reserve(3);
95 vector<BReal> 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 BReal Bernstein(int k, int n, const BReal & u)
142 return BReal(BinomialCoeff(n, k)) * Power(u, k) * Power(BReal(1.0) - u, n-k);
147 * Returns the parametric parameter at the turning point(s)
148 * In one coordinate direction
151 pair<BReal, BReal> BezierTurningPoints(const BReal & p0, const BReal & p1, const BReal & p2, const BReal & p3)
154 if (p1 == p2 && p2 == p3)
156 return pair<BReal,BReal>(0, 1);
158 BReal a = ((p1-p2)*3 + p3 - p0);
159 BReal b = (p2 - p1*2 + p0)*2;
164 return pair<BReal, BReal>(0,1);
168 return pair<BReal, BReal>(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<BReal, BReal>(0,1);
176 vector<BReal> tsols = SolveQuadratic(a, b, c);
177 if (tsols.size() == 1)
178 return pair<BReal,BReal>(tsols[0], tsols[0]);
179 else if (tsols.size() == 0)
180 return pair<BReal, BReal>(0,1);
182 return pair<BReal,BReal>(tsols[0], tsols[1]);
186 inline bool CompBRealByPtr(const BReal * a, const BReal * b)
192 * Get top most *point* on Bezier curve
194 pair<BReal,BReal> Bezier::GetTop() const
196 pair<BReal, BReal> tsols = BezierTurningPoints(y0,y1,y2,y3);
197 BReal tx0; BReal ty0;
198 BReal tx1; BReal ty1;
199 Evaluate(tx0, ty0, tsols.first);
200 Evaluate(tx1, ty1, tsols.second);
201 vector<const BReal*> v(4);
206 sort(v.begin(), v.end(), CompBRealByPtr);
207 pair<BReal,BReal> 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<BReal,BReal> Bezier::GetBottom() const
233 pair<BReal, BReal> tsols = BezierTurningPoints(y0,y1,y2,y3);
234 BReal tx0; BReal ty0;
235 BReal tx1; BReal ty1;
236 Evaluate(tx0, ty0, tsols.first);
237 Evaluate(tx1, ty1, tsols.second);
238 vector<const BReal*> v(4);
243 sort(v.begin(), v.end(), CompBRealByPtr);
244 pair<BReal,BReal> 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<BReal,BReal> Bezier::GetLeft() const
270 pair<BReal, BReal> tsols = BezierTurningPoints(x0,x1,x2,x3);
271 BReal tx0; BReal ty0;
272 BReal tx1; BReal ty1;
273 Evaluate(tx0, ty0, tsols.first);
274 Evaluate(tx1, ty1, tsols.second);
275 vector<const BReal*> v(4);
280 sort(v.begin(), v.end(), CompBRealByPtr);
281 pair<BReal,BReal> 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<BReal,BReal> Bezier::GetRight() const
308 pair<BReal, BReal> tsols = BezierTurningPoints(x0,x1,x2,x3);
309 BReal tx0; BReal ty0;
310 BReal tx1; BReal ty1;
311 Evaluate(tx0, ty0, tsols.first);
312 Evaluate(tx1, ty1, tsols.second);
313 vector<const BReal*> v(4);
318 sort(v.begin(), v.end(), CompBRealByPtr);
319 pair<BReal,BReal> 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<BReal> Bezier::SolveXParam(const BReal & x) const
343 BReal c((x1 - x0)*BReal(3));
344 BReal b((x2 - x1)*BReal(3) - c);
345 BReal a(x3 -x0 - c - b);
346 vector<BReal> 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<BReal> Bezier::SolveYParam(const BReal & y) const
359 BReal c((y1 - y0)*BReal(3));
360 BReal b((y2 - y1)*BReal(3) - c);
361 BReal a(y3 -y0 - c - b);
362 vector<BReal> 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<BReal> & 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 BRect Bezier::SolveBounds() const
387 pair<BReal, BReal> tsols = BezierTurningPoints(x0, x1, x2, x3);
389 BReal tp0; BReal tp1; BReal 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 BReal*> v(4);
401 // Not using a lambda to keep this compiling on cabellera
402 sort(v.begin(), v.end(), CompBRealByPtr);
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(), CompBRealByPtr);
422 result.h = *(v[3]) - result.y;
424 //Debug("Solved Bezier %s bounds as %s", Str().c_str(), result.Str().c_str());