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);
17 if (a == BReal(0) && b != BReal(0))
19 roots.push_back(-c/b);
22 BReal disc(b*b - BReal(4)*a*c);
27 else if (disc == BReal(0))
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 < BReal(0) && u < BReal(0)) || (l > BReal(0) && u > BReal(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*BReal(3), b*BReal(2), c));
96 //Debug("%u turning points", turns.size());
97 for (unsigned i = 1; i < turns.size(); ++i)
100 tu = std::min(turns[i],tu);
101 CubicSolveSegment(roots, a, b, c, d, tl, tu,delta);
107 CubicSolveSegment(roots, a, b, c, d, tl, tu,delta);
114 * Use dynamic programming / recursion
118 static unordered_map<int, int> dp;
119 static bool init = false;
125 auto it = dp.find(n);
128 int result = n*Factorial(n-1);
134 * Binomial coefficients
136 int BinomialCoeff(int n, int k)
138 return Factorial(n) / (Factorial(k) * Factorial(n-k));
142 * Bernstein Basis Polynomial
144 BReal Bernstein(int k, int n, const BReal & u)
146 return BReal(BinomialCoeff(n, k)) * Power(u, k) * Power(BReal(1.0) - u, n-k);
151 * Returns the parametric parameter at the turning point(s)
152 * In one coordinate direction
155 pair<BReal, BReal> BezierTurningPoints(const BReal & p0, const BReal & p1, const BReal & p2, const BReal & p3)
158 if (p1 == p2 && p2 == p3)
160 return pair<BReal,BReal>(0, 1);
162 BReal a = ((p1-p2)*BReal(3) + p3 - p0);
163 BReal b = (p2 - p1*BReal(2) + p0)*BReal(2);
168 return pair<BReal, BReal>(0,1);
170 if (t > BReal(1)) t = 1;
171 if (t < BReal(0)) t = 0;
172 return pair<BReal, BReal>(t, t);
174 //Debug("a, b, c are %f, %f, %f", Float(a), Float(b), Float(c));
175 if (b*b - a*c*BReal(4) < BReal(0))
177 //Debug("No real roots");
178 return pair<BReal, BReal>(0,1);
180 vector<BReal> tsols = SolveQuadratic(a, b, c);
181 if (tsols.size() == 1)
182 return pair<BReal,BReal>(tsols[0], tsols[0]);
183 else if (tsols.size() == 0)
184 return pair<BReal, BReal>(0,1);
186 return pair<BReal,BReal>(tsols[0], tsols[1]);
190 inline bool CompBRealByPtr(const BReal * a, const BReal * b)
196 * Get top most *point* on Bezier curve
198 pair<BReal,BReal> Bezier::GetTop() const
200 pair<BReal, BReal> tsols = BezierTurningPoints(y0,y1,y2,y3);
201 BReal tx0; BReal ty0;
202 BReal tx1; BReal ty1;
203 Evaluate(tx0, ty0, tsols.first);
204 Evaluate(tx1, ty1, tsols.second);
205 vector<const BReal*> v(4);
210 sort(v.begin(), v.end(), CompBRealByPtr);
211 pair<BReal,BReal> result;
212 result.second = *v[0];
217 else if (v[0] == &y3)
221 else if (v[0] == &ty0)
225 else if (v[0] == &ty1)
233 * Get bottom most *point* on Bezier curve
235 pair<BReal,BReal> Bezier::GetBottom() const
237 pair<BReal, BReal> tsols = BezierTurningPoints(y0,y1,y2,y3);
238 BReal tx0; BReal ty0;
239 BReal tx1; BReal ty1;
240 Evaluate(tx0, ty0, tsols.first);
241 Evaluate(tx1, ty1, tsols.second);
242 vector<const BReal*> v(4);
247 sort(v.begin(), v.end(), CompBRealByPtr);
248 pair<BReal,BReal> result;
249 result.second = *v[3];
254 else if (v[3] == &y3)
258 else if (v[3] == &ty0)
262 else if (v[3] == &ty1)
270 * Get left most *point* on Bezier curve
272 pair<BReal,BReal> Bezier::GetLeft() const
274 pair<BReal, BReal> tsols = BezierTurningPoints(x0,x1,x2,x3);
275 BReal tx0; BReal ty0;
276 BReal tx1; BReal ty1;
277 Evaluate(tx0, ty0, tsols.first);
278 Evaluate(tx1, ty1, tsols.second);
279 vector<const BReal*> v(4);
284 sort(v.begin(), v.end(), CompBRealByPtr);
285 pair<BReal,BReal> result;
286 result.first = *v[0];
291 else if (v[0] == &x3)
295 else if (v[0] == &tx0)
299 else if (v[0] == &tx1)
308 * Get left most *point* on Bezier curve
310 pair<BReal,BReal> Bezier::GetRight() const
312 pair<BReal, BReal> tsols = BezierTurningPoints(x0,x1,x2,x3);
313 BReal tx0; BReal ty0;
314 BReal tx1; BReal ty1;
315 Evaluate(tx0, ty0, tsols.first);
316 Evaluate(tx1, ty1, tsols.second);
317 vector<const BReal*> v(4);
322 sort(v.begin(), v.end(), CompBRealByPtr);
323 pair<BReal,BReal> result;
324 result.first = *v[3];
329 else if (v[3] == &x3)
333 else if (v[3] == &tx0)
337 else if (v[3] == &tx1)
344 vector<BReal> Bezier::SolveXParam(const BReal & x) const
347 BReal c((x1 - x0)*BReal(3));
348 BReal b((x2 - x1)*BReal(3) - c);
349 BReal a(x3 -x0 - c - b);
350 vector<BReal> results(SolveCubic(a, b, c, d));
351 for (unsigned i = 0; i < results.size(); ++i)
354 Evaluate(p.x, p.y, results[i]);
360 vector<BReal> Bezier::SolveYParam(const BReal & y) const
363 BReal c((y1 - y0)*BReal(3));
364 BReal b((y2 - y1)*BReal(3) - c);
365 BReal a(y3 -y0 - c - b);
366 vector<BReal> results(SolveCubic(a, b, c, d));
367 for (unsigned i = 0; i < results.size(); ++i)
370 Evaluate(p.x, p.y, results[i]);
375 vector<Vec2> Bezier::Evaluate(const vector<BReal> & u) const
377 vector<Vec2> result(u.size());
378 for (unsigned i = 0; i < u.size(); ++i)
380 Evaluate(result[i].x, result[i].y, u[i]);
386 * Get Bounds Rectangle of Bezier
388 BRect Bezier::SolveBounds() const
391 pair<BReal, BReal> tsols = BezierTurningPoints(x0, x1, x2, x3);
393 BReal tp0; BReal tp1; BReal o;
394 Evaluate(tp0, o, tsols.first);
395 Evaluate(tp1, o, tsols.second);
397 //Debug("x: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
399 vector<const BReal*> v(4);
405 // Not using a lambda to keep this compiling on cabellera
406 sort(v.begin(), v.end(), CompBRealByPtr);
409 result.w = *(v[3]) - result.x;
411 // Do the same thing for y component (wow this is a mess)
412 tsols = BezierTurningPoints(y0, y1, y2, y3);
413 Evaluate(o, tp0, tsols.first);
414 Evaluate(o, tp1, tsols.second);
417 //Debug("y: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
423 sort(v.begin(), v.end(), CompBRealByPtr);
426 result.h = *(v[3]) - result.y;
428 //Debug("Solved Bezier %s bounds as %s", Str().c_str(), result.Str().c_str());