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)
99 //if (tl > max) break;
101 CubicSolveSegment(roots, a, b, c, d, tl, tu,delta);
105 CubicSolveSegment(roots, a, b, c, d, tl, tu,delta);
111 * Use dynamic programming / recursion
115 static unordered_map<int, int> dp;
116 static bool init = false;
122 auto it = dp.find(n);
125 int result = n*Factorial(n-1);
131 * Binomial coefficients
133 int BinomialCoeff(int n, int k)
135 return Factorial(n) / (Factorial(k) * Factorial(n-k));
139 * Bernstein Basis Polynomial
141 BReal Bernstein(int k, int n, const BReal & u)
143 return BReal(BinomialCoeff(n, k)) * Power(u, k) * Power(BReal(1.0) - u, n-k);
148 * Returns the parametric parameter at the turning point(s)
149 * In one coordinate direction
152 pair<BReal, BReal> BezierTurningPoints(const BReal & p0, const BReal & p1, const BReal & p2, const BReal & p3)
155 if (p1 == p2 && p2 == p3)
157 return pair<BReal,BReal>(0, 1);
159 BReal a = ((p1-p2)*BReal(3) + p3 - p0);
160 BReal b = (p2 - p1*BReal(2) + p0)*BReal(2);
165 return pair<BReal, BReal>(0,1);
167 if (t > BReal(1)) t = 1;
168 if (t < BReal(0)) t = 0;
169 return pair<BReal, BReal>(t, t);
171 //Debug("a, b, c are %f, %f, %f", Float(a), Float(b), Float(c));
172 if (b*b - a*c*BReal(4) < BReal(0))
174 //Debug("No real roots");
175 return pair<BReal, BReal>(0,1);
177 vector<BReal> tsols = SolveQuadratic(a, b, c);
178 if (tsols.size() == 1)
179 return pair<BReal,BReal>(tsols[0], tsols[0]);
180 else if (tsols.size() == 0)
181 return pair<BReal, BReal>(0,1);
183 return pair<BReal,BReal>(tsols[0], tsols[1]);
187 inline bool CompBRealByPtr(const BReal * a, const BReal * b)
193 * Get top most *point* on Bezier curve
195 pair<BReal,BReal> Bezier::GetTop() const
197 pair<BReal, BReal> tsols = BezierTurningPoints(y0,y1,y2,y3);
198 BReal tx0; BReal ty0;
199 BReal tx1; BReal ty1;
200 Evaluate(tx0, ty0, tsols.first);
201 Evaluate(tx1, ty1, tsols.second);
202 vector<const BReal*> v(4);
207 sort(v.begin(), v.end(), CompBRealByPtr);
208 pair<BReal,BReal> result;
209 result.second = *v[0];
214 else if (v[0] == &y3)
218 else if (v[0] == &ty0)
222 else if (v[0] == &ty1)
230 * Get bottom most *point* on Bezier curve
232 pair<BReal,BReal> Bezier::GetBottom() const
234 pair<BReal, BReal> tsols = BezierTurningPoints(y0,y1,y2,y3);
235 BReal tx0; BReal ty0;
236 BReal tx1; BReal ty1;
237 Evaluate(tx0, ty0, tsols.first);
238 Evaluate(tx1, ty1, tsols.second);
239 vector<const BReal*> v(4);
244 sort(v.begin(), v.end(), CompBRealByPtr);
245 pair<BReal,BReal> result;
246 result.second = *v[3];
251 else if (v[3] == &y3)
255 else if (v[3] == &ty0)
259 else if (v[3] == &ty1)
267 * Get left most *point* on Bezier curve
269 pair<BReal,BReal> Bezier::GetLeft() const
271 pair<BReal, BReal> tsols = BezierTurningPoints(x0,x1,x2,x3);
272 BReal tx0; BReal ty0;
273 BReal tx1; BReal ty1;
274 Evaluate(tx0, ty0, tsols.first);
275 Evaluate(tx1, ty1, tsols.second);
276 vector<const BReal*> v(4);
281 sort(v.begin(), v.end(), CompBRealByPtr);
282 pair<BReal,BReal> result;
283 result.first = *v[0];
288 else if (v[0] == &x3)
292 else if (v[0] == &tx0)
296 else if (v[0] == &tx1)
305 * Get left most *point* on Bezier curve
307 pair<BReal,BReal> Bezier::GetRight() const
309 pair<BReal, BReal> tsols = BezierTurningPoints(x0,x1,x2,x3);
310 BReal tx0; BReal ty0;
311 BReal tx1; BReal ty1;
312 Evaluate(tx0, ty0, tsols.first);
313 Evaluate(tx1, ty1, tsols.second);
314 vector<const BReal*> v(4);
319 sort(v.begin(), v.end(), CompBRealByPtr);
320 pair<BReal,BReal> result;
321 result.first = *v[3];
326 else if (v[3] == &x3)
330 else if (v[3] == &tx0)
334 else if (v[3] == &tx1)
341 vector<BReal> Bezier::SolveXParam(const BReal & x) const
344 BReal c((x1 - x0)*BReal(3));
345 BReal b((x2 - x1)*BReal(3) - c);
346 BReal a(x3 -x0 - c - b);
347 vector<BReal> results(SolveCubic(a, b, c, d));
348 for (unsigned i = 0; i < results.size(); ++i)
351 Evaluate(p.x, p.y, results[i]);
357 vector<BReal> Bezier::SolveYParam(const BReal & y) const
360 BReal c((y1 - y0)*BReal(3));
361 BReal b((y2 - y1)*BReal(3) - c);
362 BReal a(y3 -y0 - c - b);
363 vector<BReal> results(SolveCubic(a, b, c, d));
364 for (unsigned i = 0; i < results.size(); ++i)
367 Evaluate(p.x, p.y, results[i]);
372 vector<Vec2> Bezier::Evaluate(const vector<BReal> & u) const
374 vector<Vec2> result(u.size());
375 for (unsigned i = 0; i < u.size(); ++i)
377 Evaluate(result[i].x, result[i].y, u[i]);
383 * Get Bounds Rectangle of Bezier
385 BRect Bezier::SolveBounds() const
388 pair<BReal, BReal> tsols = BezierTurningPoints(x0, x1, x2, x3);
390 BReal tp0; BReal tp1; BReal o;
391 Evaluate(tp0, o, tsols.first);
392 Evaluate(tp1, o, tsols.second);
394 //Debug("x: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
396 vector<const BReal*> v(4);
402 // Not using a lambda to keep this compiling on cabellera
403 sort(v.begin(), v.end(), CompBRealByPtr);
406 result.w = *(v[3]) - result.x;
408 // Do the same thing for y component (wow this is a mess)
409 tsols = BezierTurningPoints(y0, y1, y2, y3);
410 Evaluate(o, tp0, tsols.first);
411 Evaluate(o, tp1, tsols.second);
414 //Debug("y: tp0 is %f tp1 is %f", Float(tp0), Float(tp1));
420 sort(v.begin(), v.end(), CompBRealByPtr);
423 result.h = *(v[3]) - result.y;
425 //Debug("Solved Bezier %s bounds as %s", Str().c_str(), result.Str().c_str());